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PREFACE 


This  mid-term  report  includes  a  discussion  of  a  timely  and  controversial 
issue,  a  continuation  of  a  fruitful  line  of  research,  and  an  introduction  of 
new  direction  of  research,  all  oriented  toward  the  modeling  of  near-field 
records  from  underground  nuclear  explosions.  In  the  first  part,  we  compare 
the  Mueller-Murphy  and  Helmberger-Hadley  effective  source  functions,  and 
discuss  their  applicability  to  near-field  waveform  modeling.  Our  goal  is  to 
attempt  to  identify  and  reconcile  the  resolvable  differences  between  these 
source  representations  so  that  consistent  yield  scaling  relations  may  be 
obtained.  The  second  part  of  this  report  includes  a  waveform  modeling  study 
aimed  at  determining  yield  scaling  relations  for  nuclear  explosions  in  dry 
tuff  and  alluvium  at  Yucca  Flats.  This  is  a  continuation  of  the  approach  we 
have  used  for  modeling  near-field  data  from  the  Amchitka  and  Pahute  Mesa  test 
sites.  An  additional  point  of  interest,  however,  is  the  apparent 
identification  of  the  elastic  radius.  Beyond  a  scaled  slant  range  of  about 
130  ra,  elastic  waveform  modeling  predicts  observed  peak  velocities  quite  well 
but  inside  that  range,  the  observations  are  overpredicted.  In  the  final  part 
of  this  report,  we  introduce  a  method  for  the  simultaneous  inversion  of 
near-field  data  for  source  and  structure  par^uneters.  This  is  desirable  in 
order  to  obtain  quantitative  information  on  the  trade-offs  between  parameters 
and  also  to  enable  the  modeling  of  near-field  data  from  sites  for  which  no 
structure  model  previously  exists.  Preliminary  tests  are  presented  for  the 
structure  inversion  portion  of  the  method  by  inverting  synthetic  data  sets 
from  simple  models. 
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PART  1 


A  COMPARISON  OF  THE  MUELLER-HURPHY 
AND  HELMBERGER-HADLEY  SOURCE  FUNCTIONS 


by 


L.  J.  Burdick  and  J.  S.  Barker 
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INTRODUCTION 

One  of  the  most  important  goals  in  the  field  of  nuclear  seismology  is  the 
development  of  an  accurately  scaled  representation  of  the  reduced  displacement 
potential.  An  early  attempt  at  accomplishing  this  was  made  by  Haskell  (1967). 
The  most  widely  accepted  and  acknowledged  representation  was  later  presented 
by  Mueller  and  Murphy  (1971).  (An  equivalent  though  mathematically  simpler 
version  of  this  representation  was  developed  by  von  Seggern  and  Blandford 
(1972).)  Haskell  (1967)  chose  the  algebraic  form  for  his  RDP  representation 
so  that  it  was  simple  and  so  that  it  did  not  predict  discontinuities  in 
displacement,  velocity  or  acceleration.  Mueller  and  Murphy  (1971)  derived 
their  representation  by  assuming  a  form  for  the  pressure  time  history  on  the 
walls  of  a  spherical  cavity.  The  function  they  adopted  for  the  pressure 
history  is  discontinuous,  and  it  predicts  velocity  pulses  that  are 
discontinuous  and  acceleration  pulses  that  are  singular,  von  Seggern  and 
Blandford  (1972)  noted  that  observed  near-field  pulses  often  have  rise  times 
so  short  that  they  can  be  closely  approximated  by  a  discontinuous  pulse.  In  a 
more  recent  modeling  study  of  near-field  velocity  records,  Helmberger  and 
Hadley  (1981)  found  that  the  velocity  records  they  were  studying  did  have  a 
resolvable  rise  time.  To  model  this  information,  they  introduced  an  RDP 
representation  which  allows  velocity  pulses  to  have  finite  rise  times.  Their 
formalism  was  later  utilized  by  Burdick,  et  al.  (1984)  and  Hartzell,  et  al. 
(1983)  in  similar  modeling  studies  for  the  same  reasons.  It  is  unfortunate 
that  more  than  one  formalism  for  representing  the  RDP  is  in  current  use  since 
comparison  of  the  results  of  different  studies  is  made  difficult.  The  purpose’ 
of  this  investigation  is  to  find  whether  there  actually  are  any  features  of 
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the  near-field  observations  which  might  tend  to  favor  either  of  the  two 
alternate  RDP  representations  over  the  other.  A  slightly  modified 
mathematical  representation  in  which  the  Mueller  and  Murphy  ( 1971 )  source 
might  be  used  in  modeling  studies  is  explored  such  that  the  rise  time  of  the 
velocity  pulse  is  resolvably  greater  than  zero. 

The  key  results  of  the  Investigation  can  be  previewed  as  follows.  It  is 
widely  believed  that  the  most  significant  difference  between  the  two  alternate 
RDP  representations  is  that  Mueller  and  Murphy's  (1971)  has  a  high  frequency 
falloff  rate  of  f  ^  and  Helmberger  and  Hadley's  (1981)  has  a  rate  of  f~^.  We 
found  that  the  scaling  relations  in  the  two  formalisms  are  such  that  the 
spectra  are  virtually  indistinguishable  for  frequencies  from  2  and  5  Hz. 
Between  5  and  10  Hz,  they  diverge  only  by  about  a  factor  of  two.  It  has  been 
found  that  the  Muel ler-Murphy  spectra  are  definitely  too  rich  in  high 
frequen.-y  energy  in  this  band  to  satisfy  the  near-field  observations. 

However,  if  a  moderately  low  Q  value  is  postulated  for  the  near-field  crust, 
the  Muel ler-Murphy  source  can  be  made  to  agree  with  the  observations.  In 
other  words,  though  the  Mueller-Murphy  velocity  pulse  is  discontinuous,  which 
disagrees  with  the  data,  the  attenuation  corrected  pulse  is  not.  The 
corrected  Mueller-Murphy  spectra  can  be  made  to  agree  with  the  data  about  as 
well  as  the  Helmberger-Hadley  spectra  from  2  to  10  Hz.  There  are  significant 
differences  in  the  spectra  of  the  two  source  representations  in  the  frequency 
band  from  0.5  to  2.0  Hz.  This  is  the  band  which  controls  the  amplitude  of  the 

short  period  teleseismic  P  waves.  If  the  Mueller-Murphy  source  is  assumed, 

* 

the  global  average  t  value  for  events  at  NTS  comes  out  to  be  slightly  larger 
than  1.0  sec.  If  the  Helmberger-Hadley  source  is  used  it  comes  out  to  be 
about  0.7  sec.  Other  lines  of  evidence  as  to  what  this  value  should  be  might 


The  algebraic  representations  for  the  RDPs  contain  various  free 
parameters  dependent  on  yield  and  depth.  Actually,  the  functional  forms  are 
similar  enough  so  that  plots  of  the  RDPs  versus  time  will  not  appear  to  be 
significantly  different  for  comparable  values  of  the  free  parameters,  von 
Seggern  and  Blandford  (1972)  discuss  how  to  relate  the  Mueller-Murphy  free 
parameters  to  those  in  other  RDP  formalisms.  The  algebraic  expression  for  the 
Mueller-Murphy  RDP  is  given  by  the  convolution  integral 


F(t)  =  sin(bt)  exp(at)/Bb 


For  a  tuff-rhyolite  emplacement  medium 
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Y  is  yield  in  kt,  h  is  source  depth  in  m,  p  is  the  density  in  g/cm^,  Vp  is  the 
compressional  velocity  at  the  source  in  km/s,  and  A  and  p  are  Lame’s 
constants*  The  expression  for  the  Helmberger-Hadley  RDP  is 


V  (t )  -  -  exp{-kt)  P(t)  ) 

P(t)  =1  f  kt  +  -  B(kt)^ 


The  Helmberger-Hadley  scaling  relationships  for  the  same  media  (events  at 
Pahute  Mesa  below  the  water  table)  are  (Hartzell,  et  al.,  1983) 


'F  =  6.72  X  10  ” 

m 

k  -  4.04  h''  •  Vy"  •  ” 
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As  discussed  in  detail  by  Lay,  et  al.  (1984),  the  choice  of  B  =  1  at  all 
yields  was  made  arbitrarily  since  it  is  very  difficult  to  resolve  its  actual 
value  from  body  wave  data  of  any  kind.  For  events  at  the  Amchitka  test  site, 
the  scaling  relations  become  (Burdick,  et  al.,  1984) 

=  9. 5  X  10®  Y/h®  *  " 

00  ' 

k  =  4.7  h°'  Vy“‘  “ 

B  =  1 

The  most  meaningful  way  to  compare  the  source  representations  is  to  compare 
the  spectra  they  predict  in  the  0.5  to  10  Hz  band.  Graphs  for  yields  of  1000 
and  100  kt  are  shown  in  Figures  1  and  2,  respectively.  The  most  striking 
feature  that  emerges  from  this  comparison  is  the  strong  aveement  between  the 
Pahute  Mesa  spectra  between  2  and  10  Hz.  The  two  different  sets  of  scaling 
laws  were  based  on  comparable  data,  but  the  data  were  analyzed  in  very 
different  ways.  The  agreement  might  thus  be  interpreted  to  mean  that  the 
source  spectrum  is  well  resolved  in  this  band.  Note  that  particularly  at  a 
yield  of  1000  kt,  the  fact  that  the  Mueller-Murphy  source  is  decaying  as  f 
and  the  Helraberger-Hadley  as  f”^  is  not  significant.  The  differences  in  the 
scaling  of  the  corner  frequency  compensate  for  the  differences  in  the  high 
frequency  slope.  The  agreement  in  the  2  to  10  Hz  band  is  not  quite  as  close 
at  100  kt  as  1000  kt,  but  the  Hartzell,  et  al.  (1983)  data  set  contained  no 
records  from  events  this  small.  Most  of  the  records  analyzed  were  from  events 
larger  than  500  kt.  In  any  event,  the  strongest  differences  between  the 
predicted  Pahute  Mesa  source  spectra  at  either  yield  are  in  the  location  of 
the  corner  frequency  and  the  level  of  the  spectral  peak.  The  Amchitka  spectra 


Figure  2.  A  comparison  of  the  Mueller-Murphy  and  Helmberger-Hadley  yield 
scaled  RDP  spectra.  This  particular  case  is  for  100  kt. 
Helmberger-Hadley  spectra  are  shown  for  both  Pahute  Mesa  and  Amchitka 


are  shown  in  Figures  1  and  2  to  illustrate  the  point  that  the 
Helmberger-Hadley  source  was  observed  to  scale  differently  at  the  two 
different  test  sites.  An  important  point  of  controversy  is  whether  the 
characteristics  of  source  spectra  of  well  coupled  events  below  the  water  table 
vary  significantly  from  site  to  site.  The  results  of  the  Hartzell,  et  al. 
(1983)  and  the  Burdick,  et  al.  (1984)  study  seem  to  indicate  that  they  do. 

This  is  not  predicted  by  the  Mueller  and  Murphy  (1971)  correction  factors  for 
varying  emplacement  media.  No  additional  light  has  been  shed  on  this  question 
in  this  study,  but  it  is  important  to  point  out  this  significant  difference 
between  the  Mueller-Murphy  and  Helmberger-Hadley  scaling  laws. 

COMPARISONS  OF  SPECTRA  AND  WAVEFORMS 

In  this  section,  we  will  compare  the  predictions  of  the  Mueller-Murphy 
and  Helmberger-Hadley  sources  to  observed  near-field  spectra  and  waveforms 
from  well-coupled  Pahute  Mesa  events.  We  find  that  the  Helmberger-Hadley 
representation  fits  better,  but  this  is  not  suprising  since  we  are  considering 
the  data  originally  used  to  constrain  the  source  model.  On  the  other  hand, 
the  Mueller-Murphy  model  was  based  on  similar  data  from  some  of  the  same 
events. 

The  types  of  comparisons  to  be  considered  in  this  section  are  between 
completely  synthetic  spectra  and  waveforms  and  their  corresponding 
observations.  Computation  of  the  synthetics  requires  a  correction  for  the 
effects  of  wave  propagation  through  the  crust.  To  calculate  these 
corrections,  we  will  use  the  crustal  model  for  Pahute  Mesa  developed  by 
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Hartzell,  et  al.  (1983)  and  a  generalized  ray  computer  algorithm.  The  model 
of  Hartzell,  et  al.  (1983)  is  compared  to  other  crust  models  for  Pahute  in 
Figure  3.  The  model  has  no  unusual  properties  and  is  very  comparable  to  those 
found  by  H2unilton  and  Healy  (1969)  and  Carroll  (1966).  The  most  important 
point  for  this  discussion  is  that  to  P  waves  in  the  frequency  band  of 
interest,  the  crust  appears  to  have  a  smooth  positive  gradient.  This  means 
that  for  the  direct  P,  geometric  ray  theory  is  relatively  accurate.  The 
beginning  of  the  Green’s  function  of  the  crust  is  a  '  ?lta  function  arriving  at 
the  retarded  time,  and  if  the  source  pulse  initially  contains  a  discontinuity 
then  the  arriving  pulse  will  as  well.  It  is  assumed  in  the  calculations  shown 
in  this  section  that  the  Q  in  the  crust  is  high  enough  so  that  anelastic 
effects  are  not  significant.  In  the  following  section,  the  Implications  of 
this  assumption  will  be  considered  in  detail. 

It  is  generally  useful  to  consider  time  domain  Information  before 
proceeding  with  a  spectral  analysis  in  order  to  illuminate  exactly  what  the 
spectra  represent.  Also,  many  significant  spectral  characteristics  can  be 
observed  in  and  understood  from  the  time  domain  pulses.  Figures  4  and  5  show 
samples  of  near-field  P  waves  from  Pahute  Mesa.  Figure  4  shows  a  record  from 
SCOTCH  (155  let  announced)  and  BOXCAR  (1300  let  announced).  Synthetics 
predicted  by  the  Helmberger-Hadley  source  are  shown  in  dashed  line  and  those 
by  Mueller-Murphy  in  dotted  line.  The  time  window  shown  covers  about  1.5 
seconds  after  the  direct  arrival,  which  is  comparable  to  the  window  used  in 
the  spectral  analysis.  An  examination  of  the  generalized  ray  synthetics  shows 
that  the  dominant  energy  in  this  window  consists  of  direct  P  and  turning  rays 
of  the  downgoing  P  type  and  the  pP  type.  The  direct  upgoing  ray  is  about 
equal  in  amplitude  to  the  turning  ray.  Near-field  converted  PS  and  the 
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Figure  3.  P  and  S  wave  velocity  structures  for  Pahute  Mesa  (from  Hartzell 
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Figure  5.  A  comparison  of  observed  (solid  line)  and  synthetic  (dashed 
Helmberger-Hadley  or  dotted  Mueller-Murphy )  near-field  P  waves  for 
ALMENDRO  at  ranges  of  3.2  (above)  and  10  )cm  (below).  The  Mueller-Murphy 
source  is  clearly  too  high  in  frequency  content.  The  trace  amplitude  is 
indicated  on  the  right. 
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Rayleigh  wave  arrive  later.  The  roost  obvious  difference  in  the  Mueller-Murphy 
and  Helmberger-Hadley  synthetics  is  that  the  forroer  predicts  a  much  shorter 
duration  first  arrival.  It  is  rouch  too  high  in  frequency  to  roatch  the 
observations.  Figure  5  shows  a  comparison  of  data  and  synthetics  for  ALMENDRO 
(566  kt  estimated)  at  two  different  ranges.  The  fact  that  the  Mueller-Murphy 
synthetics  are  too  high  in  frequency  does  not  appear  to  depend  on  the  range. 

At  the  ranges  less  than  5  km,  the  pulses  are  relatively  simple  and  propagation 
effects  are  easily  modeled,  so  spectral  comparisons  of  data  and  synthetics 
should  yield  meaningful  results.  At  the  larger  ranges  in  Figures  4  and  5, 
multiple  reverberations  become  more  important.  They  are  more  unstable  and 
difficult  to  model  so  the  match  between  data  and  synthetics  is  poorer. 

Since  the  first  positive  peak  is  always  clear  in  the  data  and  can 
unambiguously  be  associated  with  the  direct  P  arrival,  it  is  instructive  to 
consider  its  absolute  amplitude  in  the  time  domain.  An  inspection  of  the 
trace  amplitudes  in  Figures  4  and  5  indicates  that  the  predicted  values  for 
the  Mueller-Murphy  source  are  generally  50  to  100%  higher  than  the  observed 
values,  as  should  be  expected  if  it  is  too  enriched  in  high  frequency  energy. 
The  complete  set  of  predicted  versus  observed  first  peak  amplitudes  as  well  as 
the  relative  percent  error  for  both  sources  are  given  in  Table  1.  We  define 
the  relative  percent  error  as 

Error  >•  200  ( synthetic-observed) /( synthetic+observed) . 

As  illustrated  in  Figure  6,  the  Mueller-Murphy  amplitudes  are  too  high  for  all 
events  considered.  The  event  average  of  the  relative  percent  errors  of  the 
Mueller-Murphy  amplitudes  are  indicated  by  squares  while  those  of  the 


Table  1  -  Predicted  and  Observed  Peak  Amplitudes 


Observed 

Helmberger-Hadley 

Mueller-Murphy 

Range 

Amplitude 

Amplitude 

Error 

Amplitude 

Error 

(km) 

(cm/s) 

( cm/s ) 

(%) 

(cm/s) 

(%) 

BOXCAR 

S-12* 

3.8 

82.11 

89.79 

8.9 

190.1 

79.3 

S-16* 

4.9 

60.90 

65.10 

6.7 

149.3 

84.1 

S-24 

7.3 

30.91 

29.24 

-5.5 

59.66 

63.5 

S-34 

10.4 

21.93 

11.95 

-58.9 

26.38 

18.4 

S-74 

22.5 

2.01 

3.81 

61  .8 

8.69 

124.8 

ALMENDRO 


L-03 

3.2 

67.98 

76.05 

11.2 

164.0 

82.8 

L-04 

3.2 

67.03 

76.05 

12.6 

164.0 

83.9 

ECHO 

3.4 

58.25 

59.30 

1  .8 

114.8 

65.4 

L-05 

5.1 

34.00 

47,85 

33.8 

112.2 

107.0 

L-02 

6.  1 

29.19 

28.78 

-1  .4 

53.53 

58.8 

L-06 

10.0 

7.41 

8.70 

16.0 

21.81 

98.6 

L-01 

12.6 

16.78 

5.65 

-99.2 

9.67 

-53.8 

MAST 

S-5 

3.65 

43.51 

45.37 

4.2 

112.0 

88.1 

S-6 

5.47 

13.94 

33.24 

81  .8 

77.01 

138.7 

S-7 

7.30 

8.  18 

17.72 

73.7 

36.04 

126.0 

SCOTCH 


S-3A 

4.  13 

20.05 

24.60 

20.4 

57.12 

96.1 

S-4 

6.06 

17.28 

17.32 

0.2 

35.13 

68.1 

INLET 

S-5 

1.63 

74.29 

78.38 

5.4 

167.5 

77.1 

S-6 

3.27 

28.17 

39.29 

33.0 

80.09 

95.9 

S-7 

6.53 

22.07 

19.90 

-10.3 

38.65 

54.6 

HALFBEAK 


S-5 

2.13 

91  .43 

78.93 

-7.3 

123.8 

30, 

.  1 

S-7 

2.13 

74.15 

78.93 

6.2 

123.8 

50, 

.2 

S-8 

2.28 

51 .62 

81 .51 

44.9 

1 16.6 

77, 

.2 

S-6 

3.05 

58.80 
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,6 
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Figure  6.  Average  relative  error  of  the  predicted  near-field  amplitude  of  the 
Mueller-Murphy  source  (squares)  and  the  Helmberger-Hadley  source  (stars) 
as  a  function  of  yield.  The  Mueller-Murphy  source  typically  predicts 
P-wave  amplitudes  which  are  50  to  100%  too  large  while  the 
Helmberger-Hadley  source  is  generally  accurate. 
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Helriher ger-Ha.iley  source  are  shown  as  stars.  The  only  significant  mismatch  of 
the  lattei'  source  is  for  MAST.  Hartzell,  et  al.  (1983)  suggest  that  this 
mismatch  is  due  to  the  fact  that  several  of  the  MAST  recordings  were  very 
noisy.  The  key  p)oint  of  this  brief  time  domain  analysis  is  that  it  is  clear 
from  this  comparison  t)iat  the  Mueller-Murphy  source  is  too  rich  in  high 
frequencies  to  satisfy  the  data. 

The  data  set  used  in  the  spectral  analysis  contained  24  records  from  6 
well-coupled  events  at  Pahute  Mesa.  The  data  were  originally  recorded  on 
magnetic  tape  by  Sandia  and  subsequently  machine  digitized  to  high  accuracy  at 
a  rate  of  1000  samples  per  second.  Thus,  there  is  little  question  as  to  the 
accuracy  of  the  spectra.  The  major  problem  in  the  analysis  is  in  choosing  a 
suitable  time  window.  We  have  selected  a  cosine  taper  window  of  2  seconds 
total  duration  with  a  0,25  second  taper  at  each  end.  There  are  1.5  seconds  of 
unaltered  signal.  In  each  example  to  be  shown,  an  extended  record  section 
with  the  window  placement  indicated  will  be  displayed  along  with  the  spectra 
so  the  reader  can  judge  how  appropriate  the  selection  of  the  positioning  of 
the  window  is. 

The  observed  and  synthetic  spectra  for  the  four  signals  shown  in  Figures 
4  and  5  are  shown  in  Figures  7  through  10.  The  BOXCAR  spectra  in  Figure  7 
illustrate  most  clearly  the  key  results  of  the  analysis.  At  frequencies 
between  1  and  5  Hz,  the  Mueller-Murphy  and  Helmberger-Hadley  spectra  agree 
with  each  other,  and  it  is  difficult  to  say  which  agrees  with  the  observed 
spectra  more  closely.  Between  5  and  10  Hz,  the  dotted  Mueller-Murphy  spectrum 
clearly  drifts  well  above  the  observed  while  the  Helmberger-Hadley  continues 
to  track  it.  Figure  8  compares  the  spectra  from  the  SCOTCH  record.  The  same 
general  points  are  true  though  perhaps  not  as  clearly.  This  example  was 


BOXCAR  S-34  UV  10. 40  KM 

SIGNAL  HINDOW  2.70  TO  4.70  SEC 


Figure  7.  A  comparison  of  the  observed  (solid  line)  versus  predicted 

(Mueller-Murphy  dotted  line  or  Helmberger-Hadley  dashed  line)  near-field 
P-wave  velocity  spectra  for  BOXCAR.  The  original  observed  and  synthetic 
time  series  are  shown  at  the  top  with  the  trace  amplitudes  indicated. 

The  positioning  and  length  of  the  signal  window  is  indicated  at  the  very 

p. 


SCOTCH  S-3fi  UV  4. 13  KM 

SIGNAL  HINDOW  1.00  TO  3.00  SEC 


Figure  8.  A  comparison  of  the  observed  (solid  line)  versus  predicted 

(Mueller-Murphy  dotted  line  or  Helmberger-Hadley  dashed  line)  near-field 
P-wave  velocity  spectra  for  SCOTCH.  The  original  observed  and  synthetic 
time  series  are  shown  at  the  top  with  the  TRACE  amplitudes  indicated. 

The  positioning  and  length  of  the  signal  window  is  indicated  at  the  very 


RLMENDRO  L-03  UVH  3.20  KM 

SIGNAL  HINDOW  0.80  TO  2.80  SEC 


Figure  9.  A  comparison  of  the  observed  (solid  line)  versus  predicted 

(Mueller-Murphy  dotted  line  or  Helmberger-Hadley  dashed  line)  near-field 
P-wave  velocity  spectra  for  ALMENDRC  at  3.2  km.  The  original  observed 
and  synthetic  time  series  are  shovm  at  the  top  with  the  trace  amplitudes 
indicated.  The  positioning  and  length  of  the  signal  window  is  indicated 
at  the  very  top. 


RLMENDRO  L-06  UVH  10.00  KM 

SIGNAL  WINOOW  2.50  TO  4.50  SEC 


gure  10.  A  comparison  of  the  observed  (solid  line)  versus  predicted 

(Mueller-Murphy  dotted  line  or  Helmberger-Hadley  dashed  line)  near-field 
P-wave  velocity  spectra  for  ALMENDRO  at  10  )tm.  The  original  observed  and 
synthetic  time  series  are  shown  at  the  top  with  the  trace  amplitudes 
indicated.  The  positioning  and  length  of  the  signal  window  is  indicated 
at  the  very  top. 
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chosen  to  illustrate  that  the  spectral  values  out  to  at  least  2  Hz  are  still 
fairly  sensitive  to  the  characteristics  of  the  window.  The  two  examples  from 
ALMENDRO  in  Figures  9  and  10  continue  to  show  that  the  Mueller-Murphy  spectra 
are  too  high  in  eunplitude  at  high  frequency,  and  furthermore  that  the  mismatch 
does  not  appear  to  depend  on  range.  The  spectra  at  the  more  distant  stations 
exhibit  more  complexity  and  the  choice  of  window  is  much  more  problematical. 

All  24  of  the  near-field  P-wave  spectra  are  shown  in  the  Appendix.  It  is 
not  fruitful  to  consider  them  on  a  case  by  case  basis  since  their  exact  form 
is  fairly  unstable.  It  is  necessary  to  find  some  objective  way  of  averaging 
the  results  to  isolate  the  stable  characteristics.  We  elected  to  average  the 
ratios  of  the  predicted/observed  spectra.  This  ratio  would  be  unity  at  all 
frequencies  if  the  synthetics  fit  the  observations  exactly.  The  results  of 
the  averaging  of  the  24  sets  of  spectra  are  shown  in  Figur’S  1 1 .  The 
Helraberger-Hadley  source  fits  better  than  the  Mueller-Murphy  at  all 
frequencies  greater  than  5  Hz.  The  error  in  the  Mueller-Murphy  source 
Increases  steadily  with  frequency.  In  the  1  to  5  Hz  band,  the  two  models  fit 
equally  well  but  again  the  details  of  the  spectral  ratio  in  this  range  are 
highly  dependent  on  the  choice  of  window.  In  Figure  12,  we  show  a  similar 
spectral  stack  but  in  this  instance  we  have  increased  the  window  length  from  2 
to  3  seconds.  The  instability  of  the  low-frequency  results  is  clear.  There 
is  some  change  in  the  high-frequency  results  as  well,  but  it  is  not 
significant.  The  Mueller-Murphy  spectral  stack  does  not  drift  upward  as 
strongly  and  the  Helraberger-Hadley  stack  drifts  slightly  downward.  The  reason 
for  this  is  that  the  window  has  been  extended  beyond  the  range  where  the 
synthetics  are  accurate  at  high  frequency.  The  ray  sum  has  been  truncated  so 
the  synthetic  has  no  high  frequency  energy  in  the  extended  window  while  the 
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Figure  11.  The  average  error  in  predicted/observed  spectra.  The  observed  and 
synthetic  near-field  P  waves  shown  in  the  appendix  were  windowed  and 
Fourier  transformed.  In  each  case,  the  ratio  of  predicted/observed 
spectra  was  formed  and  the  average  taken.  The  Mueller-Murphy  ratios  show 
a  significant  upward  drift  between  5  and  10  hz. 


STROKED  SPECTRAL  RATIO 


Figure  12.  The  average  error  in  predicted/observed  spectra  for  an  extended 
window  length.  This  is  the  seune  error  measure  displayed  in  Figure  11, 
but  the  seunpling  window  has  been  extended  from  2  to  3  sec. 


26 


observations  do.  An  examination  of  the  time  domain  signals  displayed  in  the 
appendix  makes  it  clear  that  extending  the  window  length  by  a  full  second 
without  adding  more  rays  is  inaccurate  at  high  frequencies  in  most  cases. 


ATTENUATION  IN  THE  CRUST 


The  preceding  discussion  has  made  it  clear  that  the  Mueller-Murphy  source 
is  too  rich  in  high  frequency  energy  to  satisfy  near-field  observations  from 
well-coupled  events  at  Pahute  Mesa  when  used  in  the  analysis  procedure 
followed  by  Helmberger  and  Hadley  (1981),  Burdick,  et  al.  (1984)  or  Hartzell, 
et  al.  (1983).  The  important  question  is,  can  the  procedure  be  modified  so 
that  the  Mueller-Murphy  scaled  sources  could  be  used?  The  key  assumption  to 
examine  is  that  the  effects  of  attenuation  are  negligible  for  near-field  P 
waves.  If  relatively  low  values  of  Q  are  assumed  for  the  crust,  attenuation 
will  obviously  remove  some  of  the  high  frequency  energy  from  the 
Mueller-Murphy  source  pulse.  To  state  the  problem  more  precisely,  can  a  Q 
distribution  of  depth,  z,  and  frequency,  f,  be  found  which  effectively 
translates  the  Mueller-Murphy  source  into  the  Helmberger-Hadley?  If  such  a 
distribution  could  be  found,  we  would  compute  from  it  an  attenuation  operator, 
A,  as  a  function  of  range,  r,  which  we  would  apply  to  our  synthetic 
seismograms  before  comparing  them  to  the  data.  The  condition  we  would  like  to 
satisfy  is 


HH(t,V)  -  A(r,Q(z,f))  *  MM(t,Y) 


where  *  represents  convolution.  A  separate  attenuation  operator  must  be 
computed  for  each  path  and  applied  to  the  appropriate  ray  before  summation. 
Since  A  is  a  smoothing  operator,  the  Helmberger-Hadley  source  is  lower 
frequency  than  the  Mueller-Murphy  and  we  have  very  little  outside  constraint 
on  Q,  it  is  very  likely  that  we  can  find  a  Q  distribution  to  at  least 
approximately  satisfy  this  condition.  The  only  point  we  will  address  in  this 
discussion  la  what  is  the  approximate  level  of  crustal  Q  required  to  make  the 
Mueller-Murphy  source  function  acceptable.  To  do  this  we  will  approximate  A 
with  a  Futterman  operator,  neglect  its  dependence  on  ray  path  and  simply  find 
an  appropriate  level  for  t*.  We  begin  by  comparing  the  spectra  of  synthetics 
for  a  1000  kt  shot  observed  at  a  range  of  10  km  as  shown  in  Figure  13  assuming 
that  the  effects  of  attenuation  are  negligible.  Based  on  the  discussion  in 
the  preceding  section,  we  can  assume  that  the  Helmberger-Hadley  synthetic 
would  fit  the  observations  whereas  the  Mueller-Murphy  synthetic  would  be  too 

high  in  time-domain  amplitude  and  too  rich  in  high  frequency.  We  found  that  a 

* 

t  value  of  0.03  sec  makes  the  spectra  of  the  two  synthetics  compatible  as 
shown  in  Figure  14.  The  travel  time  to  a  range  of  10  km  is  very  nearly  3  sec 
so  the  average  Q  of  the  crust  would  be  100.  If  the  average  Q  is  as  high  as 
300  (t*  >■  0.01  sec),  the  effect  of  the  attenuation  is  definitely  not  strong 
enough  to  make  the  Mueller-Murphy  source  acceptable.  If  we  perform  a  similar 

calculation  at  3  km,  as  shown  in  Figures  15  and  16,  where  travel  time  is  about 

* 

1  sec  (t  «  0.01  sec  for  Q  =>  100)  the  effect  of  attenuation  is  again 
definitely  not  strong  enough.  In  other  words,  the  range  dependence  of  the 
misfit  of  the  Mueller-Murphy  source  is  not  consistent  with  a  constant  average 
crustal  Q  value.  This  is  further  illustrated  in  Figure  17  where  we  plot  the 
misfit  in  peak  time  domain  amplitude  of  all  of  the  observations  of  both  the 


1000  KT  ELASTIC  VERT  10. 00  KM 

SIGNAL  WINDOW  2.50  TO  5.50  SEC 
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Figure  13.  Synthetic  spectra  for  a  Helmberger-Hadley  (solid  line)  and  a 

Mueller-Murphy  (dotted  line)  source  in  the  Pahute  Mesa  crustal  structure 
at  a  range  of  10  )cin.  Yield  Is  assumed  to  be  1000  kt.  Crustal 
attenuation  is  assumed  to  be  negligible.  The  Mueller-Murphy  spectrum  is 
significantly  higher  at  the  higher  frequencies,  and  the  amplitude  in  the 
time  domain  is  also  higher. 


1000  KT  Q«100  VERT  10.00  KM 

SIGNAL  WINDOW  2.50  TO  5.50  SEC 


Figure  14.  Synthetic  spectra  for  a  Helmberger-Hadley  (solid  line)  and  a 

Mueller-Murphy  (dotted  line)  source  in  the  Pahute  Mesa  crustal  structure 
at  a  range  of  10  )«n.  Yield  is  assumed  to  be  1000  kt.  Average  crustal  Q 
is  assumed  to  be  100.  The  spectra  and  time  domain  amplitudes  now  agree. 
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Figure  15.  Synthetic  spectra  for  a  Helniberger-Hadley  (solid  line)  and  a 

Mueller-Murphy  (dotted  line)  source  in  the  Pahute  Mesa  crustal  structure 
at  a  range  of  3  )?m.  Yield  is  assumed  to  be  1000  )ct.  Crustal  attenuation 
is  assumed  to  be  negligible.  The  Mueller-Murphy  spectrum  is 
significantly  higher  at  the  higher  frequencies,  and  the  eunplitude  in  the 
time  domain  is  also  higher. 


•.  iiu.\  a-'  1.  ^  . 


*.  fc-T. 


Relative  Error  (%) 


Figure  17.  Relative  error  of  the  predicced  versus  observed  time-doraain 

amplitudes  of  the  near-field  p  waves  in  the  data  set.  Mueller-Murphy 
values  are  shown  as  squares  and  Helmberger-Hadley  as  stars. 
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Helniberger-Hadley  and  the  Mueller-Murphy  source  as  a  function  of  range.  The 
former  fits  at  all  ranges  under  the  assumption  of  negligible  attenuation  while 
the  latter  is  about  75%  too  high  at  all  ranges.  This  means  that  if  the 
Mueller-Murphy  source  is  to  be  used,  a  Q  distribution  in  the  crust  must  be 
designed  which  causes  a  range  independent  reduction  in  frequency  content  and 
time-domain  aunplitude.  Such  a  distribution  will  probably  have  to  be  fairly 
elaborate  in  detail,  but  if  one  knew  a  priori  that  the  Mueller-Murphy  source 
was  correct,  a  suitable  distribution  could  probably  be  found  without  much 
difficulty. 

DISCUSSION 

It  is  unfortunate  that  the  results  of  this  brief  study  hinge  on  the  value 
of  Q  in  the  crust  since  it  is  such  a  poorly  known  quantity.  If  the  value  is 
in  the  range  of  50  to  100,  which  is  not  unreasonable,  then  crustal  attenuation 
has  a  major  effect  and  the  Mueller-Murphy  source  is  acceptable.  If  it  is  as 
high  as  300,  which  is  also  reasonable,  then  the  effect  is  much  smaller. 

Because  of  this  large  uncertainty,  we  have  no  clear  reason  for  favoring  one 
source  representation  over  the  other.  Additional  types  of  evidence  will  need 
to  be  considered. 

One  important  additional  constraint  on  the  accuracy  of  nuclear  event 
source  representations  comes  from  teleseismic,  short  period  P-wave  data.  The 
amplitudes  of  teleseismic  P  waves  scatter  considerably,  but  the  average  value 

for  an  event  appears  to  be  a  stable  and  meaningful  qxiantity.  Of  course,  the 

* 

average  value  of  t  in  the  earth  needs  to  be  known  to  correct  for  attenuation 


34 


along  the  propagation  path.  What  an  appropriate  average  value  for  t*  is  has 
been  a  subject  of  debate  for  some  years,  though  a  consensus  seems  to  be 
emerging.  It  is  important  to  always  associate  values  of  t*  with  a  frequency 
range  since  it  is  well  established  that  Q  does  vary  with  frequency  within  the 
body  wave  band.  At  very  long  periods,  such  as  those  observed  in  multiple  ScS 
waves  the  value  of  t*  for  P  waves  exceeds  ^  sec.  It  decreases  to  less  than  1 
sec  at  frequencies  higher  than  0.5  Hz.  In  an  attempt  to  utilize  this 
additional  constraint,  we  measured  amplitudes  of  the  short  period  P  waves  from 
the  Pahute  Mesa  events  for  which  we  had  analyzed  the  near-field  data.  We  used 
ab  amplitudes  (first  peak  to  first  trough)  since  these  are  known  to  be 
uncontaminated  by  pP.  We  translated  the  observed  amplitude  data  into  a  value 
for  m^  and  then  estimated  t*  values  on  an  event  by  event  basis.  The  results 
are  summarized  in  Table  2.  On  average,  the  Mueller-Murphy  RDP  predicts  a  t* 
value  for  teleseismic  short  period  P  slightly  in  excess  of  1  sec  and  the 
Helmberger-Hadley  predicts  a  value  of  0.73  sec.  Again  the  Mueller-Murphy  RDP 
appears  to  be  too  rich  in  high  frequency  energy.  Further  research  on 
teleseismic  and  crustal  attenuation  will  be  needed  before  a  firm  conclusion 
that  the  Mueller-Murphy  source  is  in  significant  error  will  be  warrented. 
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Table  2  -  t*  Estimates 

Event  Mueller-Murphy  Helmberger-Hadley 


t  (sec) 

t  (sec) 

BOXCAR 

1.10 

0.85 

ALMENDRO 

0.95 

0.70 

MAST 

1.00 

0.75 

HALFBEAK 

0.95 

0.70 

INLET 

1.00 

0.70 

SCOTCH 

1.10 

0.65 

average 

1.02 

0.73 
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APPENDIX 

The  following  figures  present  the  observed  and  synthetic  velocity 
waveforms  and  amplitude  spectra  for  the  24  near-field  records  analyzed  in  this 
study.  In  each  figure,  the  event  name,  station  identification,  gage  type 
(vertical  velocity  in  all  cases),  and  station  range  are  indicated  at  the  top 
of  the  plot.  The  window  applied  to  the  time  domain  waveforms  is  indicated 
above  the  waveforms.  The  window  is  2  seconds  in  total  duration  with  a  0.25 
sec  cosine  taper  applied  at  each  end,  so  the  length  of  unaltered  signal  being 
Fourier  transformed  is  1.5  sec.  The  spectra  are  direct  transforms  of  the 
velocity  waveforms,  so  are  true  velocity  amplitude  spectra.  Because  of 
windowing  effects  at  low  frequencies  and  noise  at  high  frequencies,  only  the 
portions  of  the  spectra  from  1  to  10  Hz  are  plotted.  The  observed  (DBS) 
waveforms  and  spectra  are  plotted  as  solid  lines,  the  Helmberger-Hadley 
synthetics  (H-H)  are  dashed  and  the  Mueller-Murphy  synthetics  (M-M)  are 
dotted.  The  peak  trace  amplitude  of  the  time  domain  waveforms  is  printed  to 
the  right  of  each  trace. 
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PART  II 

MODELING  NEAR-FIELD  OBSERVATIONS  FROM 
NUCLEAR  EXPLOSIONS  IN  DRY  TUFF  AND  ALLUVIUM 

by 

Terry  C.  Wallace  and  Jeffrey  S.  Barker 
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INTRODUCTION 

The  determination  of  source  parcuneters  from  any  seismic  waveform  amounts 
to  separating  source  excitation  and  propagation  phenomena.  For  underground 
nuclear  explosions,  seismic  yield  determination  is  usually  based  on  far-field 
measurements  such  as  short-period  P-wave  cunplitudes.  The  propagation  or  path 
effects  are  accounted  for  by  correcting  for  geometric  spreading  and 
attenuation.  Unfortunately,  lateral  variations  in  earth  structure  (in 
particular,  anelastic  attenuation  and  scattering)  can  introduce  considerable 
fluctuation  in  far-field  amplitudes.  Lay,  et  al.  (1984)  have  shown  that  for 
Pahute  Mesa  explosions,  there  is  an  order  of  magnitude  variation  in  the 
corrected  amplitude  of  the  first  cycle  (ab)  of  short-period  P  waves  with 
azimuth.  This  variation  in  amplitude,  which  is  due  to  pr  ipagation  effects, 
results  in  uncertainty  in  the  source  model. 

One  way  to  minimize  the  variations  in  path  effects  is  to  significantly 
reduce  the  path  lengths.  Helmberger  and  Hadley  (1981),  Stump  and  Johnson 
(1984)  and  Burdick,  et  al.  (1984)  have  shown  that  it  is  possible  to  obtain 
consistent,  detailed  explosion  source  histories  from  near-field  seismograms, 
especially  if  there  is  a  priori  information  about  the  velocity  structure. 
Near-field  records  also  have  the  advantage  that  they  allow  a  high  frequency 
representation  of  the  source  and  it  is  possible  to  separate  and  account  for 
effects  such  as  spall  and  slap  down. 

The  major  criticism  of  modeling  near-field  records  concerns  the 
assumption  that  the  material  response  is  linear.  In  the  immediate  vicinity  of 
an  explosion  the  rock  is  vaporized  and  crushed.  The  Induced  stresses  are  as 
high  as  a  megabar.  The  spherical  spreading  of  the  stress  pulse  results  in  a 
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reduction  of  the  strain  field,  and  at  some  distance  (the  elastic  radius),  the 
material  will  begin  to  behave  like  a  linear,  elastic  solid.  At  distances 
beyond  the  elastic  radius,  elastic  theory  can  be  used  to  recover  an  effective 
source  function.  Seismograms  from  inside  the  elastic  radius  will  be  grossly 
inconsistent  with  predictions  based  on  elastic  theory.  Failure  to  recognize 
the  elastic  radius  results  in  misinterpretation  of  the  source.  The  location 
of  the  elastic  radius  is  still  a  subject  of  considerable  controversy  (Trulio, 
1978;  Rodean,  1981;  Minster  and  Day,  1985),  and  can  Isest  be  determined  by 
modeling  of  very  well  Instrumented  explosions. 

In  this  section  of  this  report,  we  model  the  particle  velocity  records 
from  surface  and  borehole  gauges  for  three  nuclear  explosions  (DISCUS  THROWER, 
MUDPACK,  and  MERLIN)  detonated  in  dry  tuff  or  alluvium.  These  events  are 
unusual  in  that  they  were  designed,  in  part,  as  experiments  to  investigate 
transient  wave  propagation  (Ferret,  1971),  so  considerable  care  was  taken  with 
regard  to  geologic  criteria  in  their  siting.  The  elastic  velocity  structure 
between  the  source  and  receivers  is  known  from  boreholes  and  is  very  nearly 
planar,  with  only  minimal  distortion  due  to  faulting.  Since  the  velocity 
structure  is  well  determined,  it  is  possible  to  remove  even  complicated 
propagation  phenomena  such  as  the  interaction  between  reflected  and  refracted 
arrivals.  All  three  events  are  very  well  instrumented  and  allow  us  to 
investigate  the  wavefield  over  distance  ranges  from  a  fraction  of  the  depth  of 
burial  to  more  than  ten  times  that  depth.  We  can  see  a  clear  definition  of  a 
zone  where  elastic  theory  gives  consistent  estimates  of  the  source.  At  slant 
ranges  of  approximately  130  ra/kt^'^^  or  less,  the  amplitude  of  the  observed 
velocity  pulse  corresponding  to  the  geometric  arrival  is  smaller  than 
expected.  The  scaled  distance  of  130  m/kt^^^  corresponds  roughly  with  the 


65 


depth  of  burial,  which  suggests  that  elastic  waveform  modeling  can  be  useful 
in  determining  the  effective  source  functions  from  surface  gauges  even  within 
the  spall  zone.  Using  the  announced  yields  of  the  three  events,  we  developed 
a  scaling  relationship  between  yield  and  explosion  moment  for  events  in  dry 
tuff  and  alluvium.  The  scaling  relationship  is,  in  turn,  used  to  estimate  the 
yields  of  three  other  Yucca  Flats  events. 


NEAR-FIELD  MODELING 


The  data  used  in  this  study  are  the  radial  and  vertical  particle  velocity 
recordings  from  both  surface  and  borehole  gauges  for  the  events  DISCUS 
THROWER,  MUDPACK  and  MERLIN.  All  three  of  these  events  are  in  the 
northeastern  part  of  NTS  ( see  Figure  1  and  Table  1 )  and  were  detonated  in  dry 
tuff  or  alluvium,  which  has  a  P-wave  velocity  of  approximately  2.1  km/s.  The 
portions  of  the  seismograms  that  are  modeled  are  those  containing  the  body 
wave  arrivals.  Figure  2  shows  a  typical  data  array;  the  surface  recordings 
from  DISCUS  THROWER.  At  the  largest  distance  range  (1340  m) ,  the  fundamental 
mode  of  the  Rayleigh  wave  arrives  0.8  seconds  after  the  first  arrival.  At  the 
closest  ranges  (15  and  122  m) ,  the  direct  P  waves  are  marked  by  an  inflection 
in  the  velocity  pulse  (the  maximum  amplitude  of  the  P  waves  is  denoted  by  an 
arrow  in  Figure  2).  The  total  record  modeled  at  these  ranges  is  only  0.2 
seconds . 

The  arrival  of  the  body  waves  is  much  clearer  on  the  acceleration  data. 

As  an  example.  Figure  3  shows  a  comparison  of  the  velocity  and  acceleration 
records  for  the  event  MUDPACK.  For  the  closest  station  (45.7  m) ,  an 


Figure  1.  Location  of  the  three  events  studied  in  this  section  of  this 
report:  MERLIN  (2/16/65),  MUDPACK  (12/16/64)  and  DISCUS  THROWER 

(5/27/66) . 
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Table  1  -  Event  Parameters 


Event 

Date 

Time 

(GMT) 

Lat. 

(°N) 

Lon. 

(°W) 

Depth 

(m) 

Ann.  Yield 
(kt) 

DISCUS  THROWER 

5/27/66 

20:00 

37.178 

116.098 

337 

21 

MUDPACK 

12/16/64 

20:10 

37.178 

116.067 

155 

2.7 

MERLIN 

2/16/65 

17.30 

37.052 

116.024 

296 

10 

inflection  in  the  broad  velocity  pulse  (with  an  amplitude  of  approximately 
2  ft/sec)  represents  the  arrival  of  the  direct  P  wave.  This  arrival  actually 
represents  a  peak  in  the  acceleration  data.  The  cause  of  the  shoulder  on  the 
peak  in  acceleration  data  is  not  known,  although  it  may  be  related  to  spall 
opening,  or  to  a  non-elastic  wavefield  due  to  the  shockfront  from  the 
explosion.  At  larger  ranges,  the  shoulder  disappears.  At  the  farthest 
station  in  Figure  3  there  is  no  shoulder;  both  downgoing  and  upgoing  energy 
are  present.  In  this  study,  we  model  the  rise  time  and  ^unplitude  of  velocity 
records.  The  acceleration  records  are  used  to  confirm  the  pe^0(  of  the 
body-wave  arrival.  It  is  obvious  from  Figure  3  that  studies  using  only  peak 
velocity  measurements  will  differ  significantly  from  our  modeling  procedure  at 
small  ranges.  We  prefer  to  isolate  identifiable  arrivals  and  see  how  well 
elastic  wave  propagation  theory  can  predict  the  observations. 

The  method  we  use  to  calculate  the  P-wave  amplitude  is  generalized  ray 
theory  and  the  Cagniard-deHoop  technique  (see  Helmberger,  1968;  Langston  and 
Helmberger,  1975;  Burdick,  et  al.,  1984).  it  is  assumed  that  the  explosion 
can  be  represented  as  a  point  source  and  that  the  geologic  structure  can  be 
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4'{t)  =  ’1'  (l-e^^{l  +  kt+  ^(kt)^  +  ^(kt)  -  B(kt)  })  (2) 

where  y^represents  the  static  value  of  the  RDP  with  dimensions  of  volume,  k 
is  the  inverse  of  the  rise  time  of  the  source  function,  and  B  is  a  measure  of 
the  overshoot.  vonSeggern  and  Blandford  (1972)  modified  Haskell's  source 
representation  by  truncating  it  to  a  quadratic,  based  on  the  observed  f“^ 
spectral  decay  of  teleseismic,  short-period  body  waves.  We  have  chosen  to  use 
a  representation  of  the  source  Introduced  by  Helraberger  and  Hadley  (1981): 

'i'(t)  =4'  f  1  -  {  1  +  kt  +  i{kt)^  -  B(kt)^  )]  (3) 

The  parameters  have  the  same  meaning  as  in  Haskell's  model,  although,  except 
for  which  is  Identical  for  any  of  the  source  representations  discussed,  the 
values  of  rise  time  and  overshoot  will  vary  between  models.  Section  I  of  this 
report  provides  justification  for  our  use  of  (3)  rather  than  a  f  model: 
both  are  consistent  with  the  data  within  the  near-field  frequency  band  (2-5 
Hz)  as  long  as  the  latter  is  used  in  conjunction  with  a  depth-dependent 
crustal  attenuation  model.  For  lack  of  sufficient  a  priori  information  on 
crustal  attenuation  at  Yucca  Plats,  we  have  chosen  to  use  (3)  as  the  simpler 
model . 

There  are  three  parameters  that  must  be  modeled  in  (3);  T  ,  k  and  B.  As 

00 

discussed  by  Burdick,  et  al.  (1982)  and  Hartzell,  et  al.  (1983),  it  is  not 
possible  to  resolve  B  from  local  strong-motion  data  because  the  short  duration 
precludes  sufficient  long-period  energy.  At  the  frequencies  of  interest,  B 
trades  off  directly  with  4'^.  Since  the  main  conclusions  we  reach  in  this 
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study  do  not  depend  on  knowing  the  absolute  value  of  the  pair,  we  assume 

a  value  for  B  =  1.  Since  the  scaling  relations  which  we  use  were  derived  with 
this  constraint,  the  errors  in  the  assumption  of  B  will  not  adversely  affect 
yield  estimates.  Lay,  et  al.  (1984)  showed  that  the  amount  of  overshoot  in 
the  RDP  decreases  with  yield  (and  hence  with  increasing  depth  of  burial)  for 
the  Amchitka  teat  site.  Using  the  Lay,  et  al.  (1984)  scaling  law  gives  a 
value  of  about  2.0  for  B  for  the  events  modeled  in  this  study.  If  one  prefers 
the  higher  value  for  B,  then  4*^  scales  down  accordingly. 

The  rise  time  par<uneter,  k,  controls  the  corner  frequency.  For  the 
smallest  event  studied  (MUDPACK),  the  corner  frequency  is  higher  than  the  data 
frequency,  and  it  is  impossible  to  resolve  k.  Hartzell,  et  al.  (1983) 
developed  scaling  laws  for  k  and  B  from  strong-motion  data  from  several  Pahute 
Mesa  events: 


k 


1  1 
4.04  h^‘^  / 


(4) 


where  h  is  the  depth  of  burial  in  meters  and  Y  is  the  yield  in  kilotons.  We 
used  (5)  to  determine  k  for  MERLIN  and  DISCUS  THROWER,  found  that  the 
agreement  with  observational ly  determined  values  was  very  consistent,  so  we 
used  (5)  to  calculate  k  for  MUDPACK. 

The  static  value  of  the  RDP  was  determined  by  comparing  the  observed 
amplitude  with  the  synthetic  responses.  The  consistency  of  4*^  determined  from 
the  different  data  is  a  measure  of  the  ability  of  linear  theory  to  predict  the 
cunplitude.  Theoretically,  the  RDP  must  satisfy  the  condition  that  it  is 
invariant  with  distance.  Therefore,  assuming  the  propagation  effects  are 
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modeled  correctly,  any  systematic  variation  in  H'^with  distance  will  reflect 
non-linear  or  non- symmetric  behavior. 

THE  DATA 

DISCUS  THROWER;  The  DISCUS  THROWER  experiment  was  detonated  May  27,  1966  at 
20:00:00.0  GMT.  The  working  point  was  at  a  depth  of  337  m  in  dry  tuff  (above 
the  water  table).  Figure  4  shows  a  geologic  cross-section  along  the  line  of 
boreholes.  There  are  three  basic  units  in  the  geologic  section:  (1)  dry 
desert  alluvium  with  an  average  thickness  of  about  170  m,  (2)  a  bedded  ashflow 
tuff,  and  (3)  a  sequence  of  Paleozoic  carbonate  rocks.  Also  shown  in  Figure  4 
are  the  locations  of  various  recording  instruments.  The  length  of  the  array 
is  1.3  km,  and  the  maximum  surface  relief  is  approximately  20  m.  Except  for  a 
minor  fault  with  a  total  offset  of  only  about  21m,  there  is  no  major 
deviation  from  planar  interfaces  along  the  contact  between  units.  DISCUS 
THROWER  was  detonated  only  30  m  from  the  carbonate-tuf f  interface. 

All  of  the  available  borehole  and  surface  gauge  data  were  modeled  for 
rise  time  and  'F^.  Table  2  gives  the  velocity  structure  used  in  calculating 
the  Green's  functions  for  the  synthetics.  The  P-wave  velocities  and  densities 
are  taken  from  Perret  and  Kimball  (1971).  The  S-wave  velocities  are 
arbitrarily  set,  with  the  exception  of  the  alluvium  layer,  which  is  adjusted 
so  that  the  velocity  model  accurately  predicts  the  travel  time  of  the  Rayleigh 
wave  to  station  12S.  Minor  adjustments  were  required  in  the  layer  thicknesses 
as  compared  to  that  observed  from  the  borehole  cuttings.  For  example,  we  have 
adjusted  the  tuff-carbonate  Interface  to  be  90  m  below  the  working  point. 


Figure  4.  Surface  and  borehole  station  locations  for  DISCUS  THROWER.  The  tuff-carbonate  interface 
marks  a  major  acoustic  impedance  contrast. 
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Table  2  -  Velocity  Structure  for  DISCUS  THROWER  and  MUDPACK 


Medium 

Vp  (km/s) 

Vg  (km/s) 

Density  {q/cn?) 

Thickness  (m) 

Alluvium 

1.7 

0.9 

1.6 

168 

Tuff 

2,1 

1.3 

1 .9 

258 

Carbonate 

5.2 

2.8 

2.7 

half space 

This  is  largely  the  result  of  the  gentle  dip  of  the  interface  along  the 
profile.  The  synthetic  modeling  of  the  12S  waveforms,  as  discussed  below, 
requires  the  greater  tuff  thickness.  Since  stations  such  as  9D  and  12E  are 
very  close  to  the  interface,  their  depths  were  adjusted  oa  the  basis  of  their 
distance  from  the  interface.  As  expected,  the  velocity  model  does  an 
excellent  job  of  predicting  the  travel  time  to  all  stations. 

The  proximity  of  the  working  point  to  an  interface  with  such  a  large 
acoustic  impedance  contrast  makes  the  modeling  difficult.  Three  basic  rays 
Influence  the  P  waveforms:  (1)  the  direct  arrival,  (2)  an  arrival  that  is 
reflected  at  the  interface,  and  (3)  an  arrival  that  is  refracted  along  the 
interface.  As  an  exeunple  of  the  importance  of  these  three  arrivals.  Figure  5 
shows  the  vertical  and  radial  observations  at  station  12S  and  the  fit  of  the 
synthetic  waveforms  decomposed  into  these  three  rays.  The  peculiar  shape  of 
the  radial  waveform  is  diagnostic.  The  small  first  arrival  is  Identified  as 
the  head  wave.  The  direct  and  reflected  arrivals  must  interfere  in  such  a  way 
as  to  give  not  two  distinct  pulses,  but  rather  the  single  large  second 
arrival.  It  is  this  modeling  result  that  required  the  adjustment  to  the  tuff 
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Figure  5.  Modeling  of  surface  station  12S  for  DISCUS  THROWER.  The  observed 
records  can  be  modeled  as  the  sum  of  three  rays:  (1)  the  direct  arrival 
(2)  the  reflection  from  the  tuff-carbonate  interface  and  (3)  the  head 
wave  along  that  interface.  The  working  point  is  so  close  to  the 
Interface  that  direct  and  reflected  arrivals  interfere. 
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layer  thickness.  The  Interaction  between  reflected  and  refracted  arrivals  can 
be  seen  in  the  borehole  data  at  station  12D  (Figure  6).  At  this  distance,  the 
station  is  just  beyond  the  crossover  distance  at  which  the  refracted  wave 
becomes  the  first  arrival.  This  interference  can  be  seen  by  the  double  peak 
in  the  first  swing  on  the  radial  component. 

The  best  average  value  for  the  inverse  rise  time  parameter,  k,  for  the 
DISCUS  THROWER  data  was  determined  to  be  16  sec”^.  Table  3  gives  the  values 
of  'I'^for  the  surface  recordings.  The  average  value  is  1.22  x  10^  cm^  with  a 

Q 

standard  deviation  of  0.32  x  10  cm  .  "J*  from  the  closest  station,  3S,  which 

00 

is  only  15  m  from  the  device  borehole,  appears  to  be  Inconsistent  with  the 
other  values.  The  average  from  the  three  farther  stations  is  1.37  x  10  cm 

Q 

with  a  standard  deviation  of  0.16  x  10  cm*'. 

There  is  considerably  more  uncertainty  in  modeling  borehole  gauges  than 
those  at  the  free  surface.  Typically,  these  waveforms  have  phases  or  arrivals 
that  are  difficult  to  explain?  also,  many  of  the  gauges  are  located  near  the 
large  velocity  contrast.  For  this  reason,  static  RDP  values  from  the  borehole 
instruments  were  determined  separately  from  those  of  the  surface  data  and 
compared.  Table  4  gives  '('^for  9  different  subsurface  stations.  There  is  a 
substantial  difference  between  the  values  determined  from  the  vertical 
component  and  those  determined  from  the  radial  component.  The  average  value 
of  H*  determined  from  the  verticals  is  1.09  x  10^  cm^  with  a  standard 

OO 

o  1 

deviation  of  0.17  x  10  cm  .  This  value  is  in  fair  agreement  with  that  from 

the  surface  recordings.  The  average  value  of  4'^determined  from  the 

9  3 

subsurface  radial  components  is  0.78  x  10  cm  with  a  standard  deviation  of 
9  3 

0.52  X  10  cm  .  For  a  given  borehole,  the  values  from  the  radial  components 
below  or  above  the  interface  are  consistent,  but  they  are  not  consistent 


Figure  6.  Vertical  and  radial  borehole  seismograms  from  site  12  for  DISCUS 
THROWER.  The  tuff-carbonate  interface  is  located  at  a  depth  between 
gages  12D  and  12E. 


79 


Table  3  -  Static  RDP  values  for  surface  recordings 
from  DISCOS  THROWER 


Station 

Vertical 

(cm^) 

Radial 

(cm^) 

3S 

0.83  X 

10^ 

0.68 

X  10^ 

5S 

1.40 

M 

1.41 

fl 

9S 

1.33 

n 

1.22 

ft 

12S 

1.66 

M 

1.21 

•1 

across  the  Interface.  The  vertical  components  of  the  velocity  do  not  appear 
as  sensitive  to  the  Interface. 

Figure  7  compares  the  predicted  and  observed  amplitudes  of  the  Initial  p 

wave  as  a  function  of  slant  range.  The  value  of  H'  used  is  an  average  from 

00 

both  surface  and  subsurface  values.  The  predictions  are  quite  good,  even 
though  the  amplitudes  change  by  a  factor  of  30  and  the  distance  varies  by  a 
factor  of  10. 


MUDPACK ;  The  MUDPACK  experiment  was  detonated  on  December  16,  1964  at 
20:10:00.0  GMT.  The  working  point  was  at  a  depth  of  155  ra  and  was  very  close 
to  that  of  DISCUS  THROWER.  It  was  assumed  that  the  velocity  structures  of  the 
two  events  were  identical.  As  with  DISCUS  THROWER,  the  event  was  well 
instrumented,  with  borehole  gauges  both  above  and  below  the  carbonate 
interface. 

Table  5  gives  the  values  of  determined  by  modeling  the  P  waves.  The 
announced  yield  of  MUDPACK  was  only  2.7  kt,  so  the  rise  time  could  be  modeled 
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Table  4  -  Static  RDP  values  for  subsurface  recordings 
from  DISCUS  THROWER 


Station 

Vertical 

(cm^) 

Radial 

( cm' 

9A 

1.18  X 

lo'^ 

1.71  X 

10^ 

9b 

1.20 

n 

1.31 

If 

9D 

1.01 

M 

0.39 

If 

9E 

0.98 

n 

0.38 

If 

9F 

0.90 

n 

0.21 

II 

12B 

1.09 

n 

0.39 

If 

12D 

1.43 

n 

0.55 

H 

12E 

1.11 

n 

1.16 

11 

12F 

0.89 

n 

0.95 

n 

with  any  value  of  the  k  parameter  greater  than  20.  For  this  reason,  we 
assumed  the  scaling  relation  of  Hartzell,  et  al.  (1983)  and  calculated  a  k 
value  of  23  sec”^.  As  with  DISCUS  THROWER,  the  values  of  H*  determined  from 

CX> 

the  vertical  component  of  the  velocity  record  are  very  consistent  and  appear 
to  be  independent  of  gauge  location,  with  the  exception  of  the  very  closest 
ranges  (e.g.,  station  N1). 

The  average  value  determined  for  from  the  surface  gauges  is  9.9  x  10^ 

3  7  3 

cm  .  The  borehole  average  is  somewhat  smaller  (9.2  x  10  cm  )  but  is 

consistent  when  the  closest  stations  are  eliminated  from  the  average.  The 

radial  gauges  show  significantly  more  scatter.  The  average  value  is  7.1  x  10^ 

cm^,  but  if  station  N1  is  eliminated  the  average  becomes  8.6  x  10^  cm^. 
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Figure  7.  A  comparison 
vertical  component 
borehole  gages  are 


of  the  observed  and  predicted  amplitude  of  the 
of  velocity  for  DISCUS  THROWER.  Both  surface  and 
Included. 


Table  5  -  Static  RDP  values  from  MUDPACK 


Surface  Gages 

Slant  Range  (m) 

Vertical 

(cm^) 

1 

162 

0.89  x 

10® 

2 

180 

0.83 

n 

3 

228 

0.97 

II 

4 

302 

1.20 

n 

Borehole  Gages 

Slant  Range  (m) 

Vertical 

(cm®) 

Radial 

(cm®) 

B4-4 

303 

10® 

0.68  X 

10® 

B4-3 

195 

0.96  X 

0.89 

It 

B4-2 

179 

0.67 

II 

0.85 

•• 

B4-1 

173 

0.65 

n 

0.89 

It 

B2-3 

279 

0.78 

n 

0.99 

If 

B2-1 

261 

1.26 

n 

1 .12 

It 

N1-4 

183 

0.39 

It 

N1-3 

131 

0.17 

II 

0.27 

II 

N1-1 

102 

0.48 

n 

0.51 

n 

MERLIN;  The  MERlilN  experiment  was  detonated  on  February  16,  1965  at  17:30:00.0 
GWT.  The  working  point  was  at  a  depth  of  296  m  in  dry  tuff,  and  was  south  of 
DISCUS  THROWER  and  MUDPACK,  but  the  borehole  cuttings  indicate  that  the 
geology  is  very  similar.  MERLIN  was  instrumented  with  two  borehole  arrays, 
but  these  were  located  very  close  to  ground  zero,  so  only  the  surface  array 
data  was  used  to  determine  H*  . 

ou 

Table  6  gives  the  values  of  4*^  determined  for  MERLIN.  At  slant  ranges 

less  than  300  m,  the  values  for  appear  to  be  underpredicted.  Beyond  300  m, 

8  3 

the  average  value  for  is  3.4  x  10  cm  .  The  radial  surface  velocity  gauges 
are  again  inconsistent  with  both  the  vertical  measurements  and  with  each 
other. 


DISCUSSION 


For  all  three  of  the  events  studied,  the  near-field  records  provided  a 
consistent  measure  of  the  effective  source  function  once  a  critical  distance 
is  reached.  Inside  this  critical  distance,  the  values  we  determine  for  are 

OO 

systematically  underestimated.  Figure  8  is  an  attempt  to  compare  all  three 
events.  For  each  station  for  a  given  explosion,  the  value  of  determined  at 
that  station  is  divided  by  the  average  for  the  event.  This  is  plotted  as  a 
function  of  scaled  distance  to  remove  the  effect  of  different  yields.  Beyond 
a  distance  of  approximately  130  scaled  meters,  the  values  of  normalized 
scatter  about  1.0  by  +20%.  At  distances  less  than  130  m/kt^"^^,  the  stations' 
values  of  the  static  RDP  is  progressively  less  than  the  large  distance  mean. 

We  interpret  this  as  the  breakdown  in  elastic  wave  propagation  theory  in 
predicting  the  ground  velocity  very  near  an  explosion.  On  the  basis  of 
non-linearity  or  shockwave  fronts,  the  underprediction  of  is 
counter-intuitive;  inside  the  elastic  radius  we  expect  the  velocity 
amplitudes  to  be  much  greater  than  predicted.  This  is  most  likely  an  artifact 
of  the  way  we  are  measuring  the  peak  velocity.  It  is  in  this  region  that  the 
acceleration  pulses  have  pronounced  shoulders,  so  the  energy  in  the  pulse 
appears  to  be  partitioned  into  at  least  two  wavefields,  only  one  of  which  is 
behaving  elastically.  If  we  look  only  at  the  peak  velocity  and  not  try  to 
identify  the  geometric  arrival,  the  values  of  tend  to  be  overpredicted. 
Although  we  do  not  have  a  model  for  the  acceleration  pulse  shoulder,  its 
presence  can  be  used  to  determine  where  elastic  waveform  modeling  may  be  used 
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Figure  8.  A  comparison  of  the  normalized determined  at 
station  divided  by  the  average  value  for  that  event)  for  the 
studied.  Beyond  a  scaled  distance  of  about  130  m/kt  ,  the 
provides  a  very  consistent  measure  of 
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Figure  9.  A  comparison  of  the  ratio  of  observed  to  synthetic  peak  velocity 
values  for  ten  Pahute  Mesa  events  modeled  by  Barker,  et  al.  (1985)  and 
Hartzell,  et  al.  (  1983).  This  measure  is  equivalent  to  the  normalized  "f, 
of  Figure  8.  For  Pahute  Mesa  events,  consistent  amplitudes  are  obtained 
beyond  a  scaled  range  of  about  175  m/kt 
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stations  have  a  horizontal  angle  of  incidence,  and  station  5S  has  a  nearly 
vertical  (up)  travel  path.  If  the  explosion  was  not  symmetric,  one  would 
expect  to  see  variations  in  the  estimate  for  H'  •  The  large  data  set  here 
argues  against  a  source  asymmetry  explanation  for  tectonic  release. 

For  all  three  events,  the  radial  P-wave  amplitudes  are  much  more 

difficult  to  predict  than  their  vertical  component  counterparts.  At  large 

ranges,  the  surface  radial  values  for  4*^  are  in  good  agreement  with  the 

vertical  estimates,  but  at  the  borehole  stations,  they  are  generally  smaller 

than  expected,  and  are  very  sensitive  to  the  location  of  the  gauge  relative  to 

the  tuff-carbonate  interface.  There  are  three  basic  hypotheses  to  explain  the 

inconsistencies:  (1)  the  horizontal  motions  are  more  influenced  by 

incorrectly  modeled  structural  pareimeters  than  the  vertical  motions,  (2)  the 

horizontal  motions  are  more  influenced  by  non-linear  behavior,  or  (3)  somehow 

the  radial  motion  is  more  difficult  to  measure  (i.e.,  poor  coupling  in  the 

borehole)  than  the  vertical  motion.  It  is  beyond  the  scope  of  this  study  to 

differentiate  between  these  possibilities,  the  variation  in  radial  amplitudes 

is  important  in  light  of  other  techniques  for  calculating  the  static  level  of 

the  RDP.  One  of  the  most  common  techniques  is  to  integrate  the  radial 

displacement  as  recorded  at  a  laorehole  gauge  at  approximately  shot  level. 

Perret  and  Kimball  (1971)  determined  'Fusing  this  technique  for  several  of  the 

borehole  stations  we  modeled  in  this  study  for  DISCUS  THROWER.  The  value  they 

9  3 

obtain  for  the  RDP  Is4.2x10  cm,  or  a  factor  of  4  larger  than  our  result. 
Figure  10  includes  a  comparison  of  the  scatter  for  the  two  methods  of 
determining  4*^.  The  integration  of  the  radial  displacement  has  scatter  of  a 
factor  of  5.0  about  the  mean,  while  using  the  elastic  modeling  of  the  vertical 
component  of  the  P  waves  has  a  scatter  of  0.2. 
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Figure  10.  A  comparison  of  methods  for  determining  4*  .  The  top  graph  shows 
the  scatter  in  as  determined  by  integration  of  the  radial  components 
of  displacement  near  shot  level.  The  lower  graph  gives  the  scatter  in  'I' 
determined  by  generalized  ray  modeling. 
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In  order  to  estimate  the  yields  of  underground  explosions  based  on  the 
effective  source  functions  as  determined  by  near-field  modeling,  we  must  have 
a  series  of  scaling  relationships  that  are  applicable  to  the  material 
properties  of  the  source  medium.  Hartzell,  et  al.  (1983)  developed  the 
following  scaling  relation  for  Pahute  Mesa  events; 


4*  *  6.72  X  10®  Y  / 

00  ' 


(5) 


where  h  is  the  depth  of  burial  in  meters.  Further,  on  the  basis  of  announced 
yields,  it  is  possible  to  write  an  empirical  relationship  for  h  and  yield: 


h  =»  0.120  Y^/^. 


(6) 


Combining  (5)  and  (6)  leads  to 

4*  -  9.6  X  10^  Y^*^. 

00 


(7) 


We  have  chosen  to  approximate  this  with  a  simple  linear  relationship  between 

yield  and  4'  ; 

00 

4*  =  10.2  X  lo’’'  Y.  (8) 

We  can  derive  the  same  sort  of  function  for  the  Yucca  Flats  data  modeled  in 
this  study  using  the  announced  yields.  This  gives 

4^  -  4.0  X  lo"^  Y. 


(9) 


This  implies  that  the  coupling  associated  with  the  dry  tuff  of  alluvium  is  2.5 
times  less  than  that  of  the  wet  tuff  of  Pahute  Mesa.  This  is  consistent  with 
the  analysis  of  Perret  and  Kimball  (1971),  who  found  that  DISCUS  THROWER  had  a 
very  low  seismic  efficiency;  only  0.25%  of  the  energy  released  in  the 
explosion  reaches  the  region  of  elastic  response  as  seismic  energy.  In 
contrast,  explosions  detonated  in  hard  rock  typically  have  an  efficiency  of 
about  one  percent. 

We  can  use  (9)  to  predict  the  yields  of  other  events  detonated  in  dry 
tuff  or  alluvium.  Perret  and  Bass  (1975)  summarize  the  peak  velocity 
measurements  for  several  Yucca  Flats  events.  We  took  the  velocity  amplitudes 
and  station  parameters,  along  with  the  velocity  model  used  for  DISCUS  THROWER 
and  determined  values  of  'i'  for  PACKARD,  FISHER  and  VULCAN.  In  Figure  11  the 
values  of  are  plotted  against  announced  yield  for  these  events.  In  all 
three  cases,  the  scaling  relationship  did  a  very  good  job  of  predicting  the 
yield.  In  fact,  the  worst  event  predicted  is  DISCUS  THROWER,  which  appears  to 
fall  between  the  curves  for  wet  and  dry  tuff.  This  may  be  the  result  of 
DISCUS  THROWER  being  detonated  so  close  to  the  carbonate  Interface. 

CONCLUSIONS 

It  is  jxDssible  to  use  the  rise  time  and  eimplitude  of  the  near-field  P 

waves  to  determine  consistent  static  levels  of  for  an  effective  source 

00 

function  described  by  Helmberger  and  Hadley  (1981).  For  the  three  events 
studied  (Table  7  summarizes  the  effective  source  functions  obtained),  the 
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Figure  11.  A  linearized  scaling  law  relating  to  yield  for  events  detonated 
in  dry  tuff  or  alluvium.  Also  shewn  is  an  .approximation  to  the  relation 
developed  by  Hartzell,  et  al.  (1983)  for  Pahu  e  Mesa  events. 
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Event 

K 

(s“'') 

B 

4* 

fCB?) 

Mo 

(dyne-cm) 

Depth 

(m) 

Ann.  Yield 
(kt) 

DISCUS  THROWER 

16 

* 

1 

1.22  X  10® 

1.02  X  10^® 

337 

21 

MUDPACK 

23* 

* 

1 

9.90  X  10^ 

8.30  X  10^® 

155 

2.7 

MERLIN 

20* 

* 

1 

3.37  X  10® 

2.82  X  10^® 

296 

10 

assumed  values. 


vertical  component  of  velocity  can  be  modeled  using  generalized  rays  for 
scaled  slant  ranges  greater  than  130  mAt^'^^.  Since  this  roughly  corresponds 
to  the  depth  of  burial,  this  implies  that  source  strength  may  be  determined 
from  surface  gauges  within  the  spall  zone  without  having  to  resort  to 
extensive  modeling  of  non-linearity  and  spall  or  slap-down. 

The  scaling  relationship  which  was  derived  relating  ’1’^  to  yield  implies 
that  coupling  is  2.5  times  smaller  for  dry  tuff  and  alluvium  shots  than  for 
the  wet  tuff  events  at  Pahute  Mesa.  The  fact  that  once  the  region  was 
calibrated  by  an  event  where  the  velocity  structure  is  known  in  detail  (such 
as  DISCUS  THROWER)  a  scaling  relationship  can  be  used  to  predict  the  yield  of 
other  events  in  the  region  can  be  of  fundamental  importance  in  threshold 


monitoring 
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PART  Ill 


SIMULTAHEOUS  INVERSlOH  OF  NEAR-FIELD  DATA 
FOR  SOURCE  AND  STRUCTURE  PARAMETERS  - 
A  PRELIMINARY  ASSESSMENT 

by 

Jeffrey  S.  Barker  and  L.  J.  Burdick 


INTRODUCTION 


The  determination  of  source  parameters  for  underground  nuclear  explosions 
through  waveform  modeling  of  near-fleld  recordings  is  dependent  upon  the 
velocity  structure  model  assumed  between  the  source  and  receivers.  While 
Independent  information  may  exist  for  some  aspects  of  the  structure,  such  as 
the  velocity  at  the  surface  and  the  depths  to  the  water  table  and  basement, 
the  detailed  models  are  generally  the  result  of  laborious  trial-and-error 
modeling  of  observed  waveforms  (e.g.,  Helmberger  and  Hadley,  1981;  Hartzell, 
et  al.,  1983;  Burdick,  et  al.,  1984).  Discretizations  of  velocity  gradients 
into  plane-layered  structure  models  must  be  fine  enough  to  avoid  introducing 
high  frequency  errors  into  the  source  function  determinations.  For  signals 
with  frequencies  up  to  5  Hz  and  compressional  velocities  near  3  km/s,  layer 
thickness  should  be  less  than  about  0.6  km.  This  is  particularly  important 
near  the  source  depth  due  to  the  Interference  of  upgoing  and  diving  energy 
over  near-fleld  ranges.  Thus  an  adequate  structure  model  for  near-field 
waveform  modeling  might  consist  of  more  than  ten  layers  in  the  top  10  km  of 
the  crust.  There  is  obviously  a  premium  to  be  paid  before  undertaking  the 
nK>deling  of  events  at  different  test  sites.  Further,  the  structure  model 
obtained  is,  to  a  certain  extent,  dependent  upon  the  source  function  used 
during  the  modeling  procedure.  The  result  is  clearly  non-unique  and  the 
extent  of  the  trade-offs  cannot  easily  be  quantified. 

One  solution  to  this  problem  is  a  simultaneous  inversion  for  source  and 
structure  parameters.  The  complexity  of  wave  interaction  in  the  near-fleld 
makes  such  an  inversion  an  approximate  and  somewhat  costly  endeavor.  So, 
while  the  models  obtained  by  such  an  inversion  may  not  be  intrinsically 


"better"  than  their  trial-and-error  counterparts,  a  suitable  choice  of 
inversion  toethod  will  automatically  provide  Invaluable  information  concerning 
parameter  and  data  resolution.  Not  only  can  the  source  vs.  structure 
trade-offs  be  quantified  in  this  way,  but  also  the  trade-offs  within  the 
source  and  structure  parameterizatlons  (such  as  the  B-4'^  trade-off  in  the  RDP) 
may  be  determined.  We  also  obtain  an  indication  of  which  aspects  of  the 
observed  data  are  most  important  in  constraining  the  solution. 

This  section  of  this  report  includes  a  discussion  of  the  development  of 
such  an  Inversion,  as  well  as  preliminary  tests  of  the  structure  model 
determination  aspect  using  synthetic  data  sets.  Application  of  the  method  to 
near-field  data  from  Pahute  Mesa  and  Yucca  Flats  events  will  be  Included  in 
the  Final  Technical  Report. 

BACKGROUND 

The  vertical  displacement  due  to  an  explosion  source  may  be  written  using 
first  order  asymptotic  Generalized  Ray  theory  as  (see  for  exeunple,  Helmberger 
and  Harkrider,  1978) 

rays 

w(r,t)  =  *  I  ^^(P)  ) 

where  *  denotes  convolution,  'l”(t)  is  the  time  derivative  of  the  reduced 
displacement  potential  (RDP),  and 
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^Ap)  =  i? )c 

a 

which  is  evaluated  along  the  Cagniard  contour,  p  is  the  ray  parameter,  n  is 

a 

the  vertical  wave  slowness  ( -  p^  ),  receiver  function 

defined  by  Helmberger  (1974),  and  nj^(p)  is  the  product  of  reflection  and 
transmission  coefficients  along  the  path  for  each  ray. 

For  the  inversion,  we  need  to  compute  the  partial  derivatives  of  ( 1 )  with 
respect  to  the  individual  parameters  of  the  source  and  velocity  structure. 

The  partial  with  respect  to  structure  parameters,  say  the  velocity  of  the  jth 
layer,  Vj,  reduces  to  the  evaluation  of 

9r^(p)  3R^(p)  9p 

9v ,  9p  9v. 

3  3 

which  consists  of  evaluating  the  partials  of  the  various  terms  in  (2)  as  a 
function  of  the  complex  ray  parameter  along  the  Cagniard  contour,  while  the 
contour  itself  varies  with  changes  in  the  velocity  parameter.  This  evaluation 
must  be  made  for  each  ray,  then  summed  to  obtain  the  partial  derivative. 

Previous  structure  inversions  based  on  waveform  modeling  have  invoked  a 
variety  of  approximations  in  order  to  simplify  (3).  Mellman  (1980)  utilized  a 
modified  first  motion  approximation,  such  that  rapidly  varying  terms  in  (2) 
are  approximated  by  evaluating  them  only  at  the  geometric  ray  parameter.  Shaw 
(1983)  and  Chapman  and  Orcutt  (1985)  used  the  WKBJ  approximation  to  invert 
marine  reflection/refraction  data,  and  Brown  (1982)  employed  a  similar 
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approach  to  obtain  the  ocean  sound  speed  profile.  Since  the  WKBJ 
approximation  becomes  invalid  where  large  velocity  gradients  exist,  or  where 
discontinuities  are  encountered  near  the  turning  depth  of  a  ray.  Given  (1984) 
developed  a  hybrid  WKBJ-Generalized  Ray  approach  to  obtain  upper  mantle 
structure.  In  Inverting  long-period,  regional  waveforms  for  average 
crustal  structure,  Wallace  (1983)  generated  average  Generalized  Ray  Green's 
functions,  which  were  "stretched"  or  "squeezed"  to  simulate  perturbations  in 
the  structure  model. 

In  forward  modeling  of  near-field  waveforms,  Burdick,  et  al.  (1982)  and 
Barker,  et  al.  (1985)  have  shown  that  the  observed  waveforms  result  from  the 
Interference  of  upgoing  and  diving  rays,  as  well  as  surface  reflections.  This 
Interference  is  particularly  sensitive  to  ti a  velocity  gradient  near  the 
source  depth.  Each  synthetic  seismogram  consists  of  a  r  'perposition  of  rays 
from  a  variety  of  take-off  angles  that  sample  different  portions  of  the 
structure.  An  algorithm  designed  to  determine  if  and  when  the  modified  first 
motion,  WKBJ  and  other  approximations  are  appropriate  for  a  particular  ray,  to 
compute  the  analytic  partial  derivatives,  and  sum  over  rays  would  seem  to  be 
of  questionable  efficiency  and  accuracy.  Instead,  numerical  partlals  may  be 
computed  by  evaluating  the  terms  in  (3)  along  the  complex  ray  parameter  for 
each  ray.  This  approach  is  obviously  nonlinear,  and  care  must  be  taken  to 
limit  the  model  perturbations  from  which  partials  are  computed  and  to  damp  the 
resulting  par^uneter  changes.  For  simpler  problems,  the  computation  of 
numerical  partial  derivatives  would  be  terribly  inefficient,  but  given  the 
complexity  of  the  near-field  waveform  modeling  problem,  the  accuracy  of  the 
method  will  more  than  counterbalance  the  computational  cost. 

A  number  of  techniques  have  been  used  to  obtain  source  parameters  by 


inversion  of  body  waveforms.  Burdick  and  Mellman  (1978)  inverted  teleseismic 
body  waves  to  obtain  the  orientation  and  time  history  of  the  Borrego  Mt. 
earthquake.  Several  investigators  have  inverted  teleseismic  body  waves  to 
obtain  the  seismic  moment  tensor.  A  few  examples  include  Stump  and  Johnson 
(1977),  Langston  (1981),  Dziewonski,  et  al.  (1980)  and  Nabelek  (1984).  To  our 
knowledge,  the  only  inversion  of  near-field  data  for  explosion  source 
par^uneters  was  performed  by  Stump  and  Johnson  (1984).  In  that  study,  they 
solved  for  the  six  Independent  elements  of  the  moment  tensor  and  interpreted 
the  relative  amplitudes  and  time  histories  of  the  isotropic  part  and 
deviatoric  part  in  terms  of  explosion  and  tectonic  release  sources, 
respectively.  Because  the  moment  tensor  was  free  to  take  on  any  value,  and 
because  the  velocity  structure  model  was  fixed,  it  is  entirely  possible  that 
errors  in  the  structure  model  and  noise  in  the  data  could  map  into  the 
non-isotropic  moment  tensor.  In  addition,  since  no  constraint  was  placed  on 
the  form  of  the  isotropic  time  history,  only  the  low-frequency  level  could  be 
interpreted  in  terms  of  the  parameters  of  an  explosion  RDP. 

Because  of  our  assertion  that  errors  in  the  velocity  structure  adversely 
affect  explosion  source  parameter  estimates,  we  are  interested  in  approaching 
these  problems  before  tackling  source  complexity.  Therefore,  our  preliminary 
inversions  will  constrain  the  explosion  source  to  be  isotropic  and  we  will 
constrain  the  source  time  history  to  have  the  form  of  the  time  derivative  of 
the  RDP.  Various  forms  of  the  RDP  will  be  considered.  After  evaluating  these 
results,  the  source  parameterization  could  be  generalized  to  allow  the  source 
time  history  to  take  on  an  arbitrary  shape.  At  a  later  stage,  a  deviatoric 
moment  tensor  source  representing  tectonic  release  could  be  added  to  the 
isotropic  source,  with  the  latter  constrained  to  the  known  depth  and  origin 
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time  of  the  explosion.  Inverting  for  a  single,  six-component  moment  tensor, 
however,  constrains  the  non-isotropic  portion  of  the  source  to  occur  at  the 
same  location  as  the  explosion.  While  this  may  turn  out  to  be  the  case, 
near-field  observations  should  allow  resolution  of  spatial  separations  in  the 
two  sources,  so  such  a  constraint  may  be  unreasonable. 

INVERSION  METHOD 

The  layered  structure  models  required  for  Generalized  Ray  synthetics  may 
be  considered  discretizations  of  one  or  more  velocity  gradients  separated  by 
first-  or  second-order  discontinuities.  In  order  to  minimize  the  number  of 
free  parameters  in  the  inversion,  we  will  utilize  this  fact  and  solve  not  for 
the  velocities  and  thicknesses  of  Individual  layers,  but  for  the  velocity 
gradient  and  the  velocity  at  the  top  of  the  gradient.  The  discretized  model 
will  follow  this  gradient,  with  layer  thicknesses  determined  by  the  wavelength 
of  the  predominant  period  in  the  data.  Shear  velocity  and  density  may  be  free 
parameters,  or  may  be  related  to  the  compresslonal  velocity  by  predefined 
constants.  In  general,  the  free  parameters  are  the  gradient  within  each 
layer,  the  velocity  at  the  top  of  the  layer,  and  the  depth  to  the  top  of  the 
layer.  Thus,  for  an  n-layered  structure  in  which  shear  velocity  and  density 
are  not  free  parameters,  the  number  of  structure  parameters  in  the  Inversion 
is  3n-1.  The  number  of  gradients  to  be  considered  will  be  fixed  for  each 
inversion,  so  the  optimal  number  will  have  to  be  determined  by  test  Inversions 
proceeding  from  simpler  to  more  complex  models.  Constraints  such  as  the  depth 
to  the  water  table,  or  that  a  particular  interface  may  have  only  a 
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second-order  discontinuity,  -reduce  the  number  of  free  parameters  while 
increasing  non-linearity  of  the  inversion. 

Since  the  discretized  velocity  model  will  vary  from  iteration  to 
iteration  as  well  as  with  the  model  perturbations  required  to  compute  the 
numerical  partials,  the  ray  set  from  which  the  synthetics  are  computed  will 
also  vary.  Barker,  et  al.  (1985)  found  that  the  first  one  or  two  cycles  of 
vertical  velocity  observations  at  near-field  ranges  were  well  modeled  by 
including  the  upgoing  direct  ray,  a  sum  of  diving  rays  that  constitute 
turning,  reflecting  and  critically  refracting  energy,  and  a  sum  of  rays  that 
depart  upward,  reflect  from  the  free  surface  and  follow  diving  ray  paths.  If 
radial  observations  are  to  be  modeled,  or  if  later  portions  of  the  waveform 
are  considered,  P-to-S  conversions  at  the  free  surface  and  at  discontinuities 
should  be  included.  An  algorithm  to  generate  this  ray  set  for  each  new 
layered  structure  is  simple  to  implement,  and  computational  time  is  saved  by 
culling  rays  which  arrive  outside  the  inversion  time  window,  or  whose 
amplitude  will  be  below  some  predefined  cutoff  level. 

The  source  parameterization  follows  the  form  of  the  RDP  proposed  by 
Haskell  (1967): 

H'  (t)  -  4'_^(1  -  [1  +  Kt  +  (Kt)^  +  (Kt)^  -  B(Kt)^  ]  )  (4) 

where  V ^represents  the  static  value,  K  represents  the  rise  time,  and  B 
represents  the  overshoot.  Helmberger  and  Hadley  (1981)  truncated  the 
polynomial  in  (4)  to  cubic,  while  vonSeggern  and  Blandford  (1972)  suggest  that 
it  should  be  quadratic.  The  effective  source  function  for  any  of  these  forms 
may  be  obtained  by  taking  the  first  (displacement),  second  (velocity)  or  third 
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(acceleration)  time  derivatives  of  the  RDP.  As  may  be  seen  from  (4),  the  RDP 
and  its  time  derivatives  are  linear  in  and  B,  but  not  in  K.  Analytic 
partial  derivatives  may  be  computed  for  and  B,  but  we  will  rely  on 
numerical  partial  derivatives  for  K.  The  computation  of  the  effective  source 
functions  is  quite  fast,  and  the  difference  may  be  computed  before  the 
convolution  in  (1).  The  first  section  of  this  report  indicated  that  and  B 
trade  off  nearly  totally  in  modeling  near-field  velocity  records.  While  this 
trade-off  may  be  quantified  by  allowing  both  parameters  to  vary  in  the 
inversion,  useful  source  estimates  may  require  that  B  be  constrained. 

The  linearized  least-squares  Inverse  may  be  obtc.ined  by  solving 

Ac-=A-Ap* 

in  which  Ac'  is  the  residual  vector,  containing  the  differences  between 
observed  and  synthetic  data.  A*  is  the  matrix  of  partial  derivatives,  andAp’ 
is  the  desired  vector  of  parameter  changes.  Perhaps  the  simplest 
determination  of  the  residual  vector  is  a  polnt-by-polnt  difference  of 
observed  and  synthetic  seismogreuns  within  the  time  window.  This  implies, 
however,  that  neighboring  points  in  a  seismogram  represent  independent 
observations,  which  is  generally  not  the  case.  On  the  other  hand,  since  each 
seismogram  reflects  a  lagged  sum  of  rays  which  sample  different  portions  of 
the  structure,  the  resolution  of  parameters  varies  through  the  time  window, 
and  a  single  observation  for  each  seismogram  underestimates  the  information 
contained  in  the  waveform. 

Burdick  and  Mellman  (1978)  utilized  an  error,  or  residual,  function  based 
on  the  normalized  cross-correlation  coefficient  between  observed  and  synthetic 


waveforms.  That  is 


=  1  -  max(Oj^  cc  s^) 


(6) 


where 


o.  (t) 
1 


o.  (t) 
1 


(Jo^(t)  dt)'- 


\/1 


(7) 


and 


s.  (t) 
1 


s.  (t) 
1 


(8) 


o^(t)  is  the  ith  observed  seismogrcun,  s^(t)  is  the  ith  synthetic  seismogram 
and  cc  denotes  cross-correlation.  This  error  function  has  become  popular 
because  it  is  sensitive  to  the  shapes  of  the  waveforms,  while  the 
normalization  makes  it  insensitive  to  absolute  amplitude,  and  measuring  the 
maximum  of  the  cross-correlation  makes  it  insensitive  to  absolute  time.  In 
the  near-field  problem,  however,  absolute  amplitude  and  time  are  well 
calibrated,  and  must  be  fit  by  the  resultant  model.  Therefore,  in  addition  to 
(6),  our  residual  will  include  for  each  seismogram  the  relative  difference  in 


the  normalization  factors. 
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(/o^(t)  dt)  -  (/s^(t)  dt) 

(/oJ(t)  dt)"/^ 


(9) 


which  is  a  measure  of  absolute  amplitude  residual,  and  the  time  lag  to  the 
maximum  of  the  cross-correlation,  tj^,  which  is  a  measure  of  absolute  time 
residual.  The  residual  vector  is  now  defined  as 

Ac-  -  [e^,  a^,  t,,  e2,  a2,  t2,  e^,  a„,  t^l^ 

for  m  seismograms  and  the  objective  function  to  be  minimized  by  the  inverse  is 

r  =  I  +  a^2  +  t^^^,  ^ 
m 

With  this  choice  for  the  residual  function,  the  numerical  partial  derivatives 
become,  rather  than  (3), 

9e  (e  (p  +Ap  )  -  e  (p  -Ap  )) 

_i  =  — i — 3 - 3 - i — 3 - 3 —  (12) 

9Pj  2  APj 

where  APj  is  the  perturbation  to  the  starting  model  parameter,  p^.  The 
partials  3a^/9pj  and  9tj^/9pj  have  similar  form. 

Now  that  we  have  defined  the  form  of  the  residual  and  pareuneter  change 
vectors,  we  may  redefine  the  matrices  in  (5)  as 


.  •  k  r  . 


in  which  s  and  w  are  the  covariance  matrices  of  the  observations  and 
parameters,  respectively.  The  basic  purpose  of  these  weighting  factors  is  to 
counter  the  dimensionality  of  the  residuals  and  the  pareuneterization  of  source 
and  structure.  For  simplicity,  we  assume  that  errors  are  iincorrelated.  If  an 
estimate  of  the  variance  is  known  {e.g.,  timing  errors),  that  value  is  used  as 
the  weighting  factor.  Otherwise,  weights  are  assigned  based  on  an  estimate  of 
the  dimensions.  For  excunple,  waveform  residuals  will  have  values  between  0.0 
and  1.0,  while  absolute  time  residuals  are  more  likely  to  be  less  than  O.ls. 
Similarly,  a  reasonable  variation  in  the  velocity  at  the  top  of  a  layer  might 
be  about  1  km/s,  while  the  variation  in  gradient  might  be  about  0.1  km/s/km. 
The  weighting  factors  may  also  be  modified  to  tailor  the  sensitivity  of  the 
inversion.  For  instance,  we  may  wish  to  preferentially  fit  absolute  travel 
time  at  the  expense  of  slight  waveform  mismatches.  In  this  case,  we  would 
decrease  the  time  residual  variance  relative  to  that  of  the  waveform  residual. 
Following  Jackson  (1979),  a  priori  constraints  may  be  placed  on  the  parMeter 
changes  so  that,  for  ex^unple,  the  velocity  at  the  top  of  a  layer  cannot  be 
negative,  or  so  that  source  pareuneters  must  fall  within  a  range  of  reasonable 
values. 

The  weighted  version  of  (5)  may  be  solved  by  a  linearized  generalized 
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inverse  (Wiggins,  1972)  such  that 


Ap*A  Ac*V/r'u'^Ac 


in  which  A  is  the  diagonal  matrix  of  non-zero  singular  values  of  A,  V  is  the 
corresponding  matrix  of  eigenvectors  which  spans  the  parouneter  space,  and  U  is 
the  matrix  of  eigenvectors  that  spans  the  residual  space.  Since  the  problem 
is  non-linear,  the  inversion  is  iterative  and  minimizes  both  the  least-squares 
error,  |a  Ap  -  Acp,  and  the  length  of  the  par^eter  changes,  |Ap|^.  Since 
small  singular  values  result  In  large  parzuneter  changes,  an  a  priori  cut-off 
is  specified,  below  which  the  singular  value  is  assumed  to  be  zero,  and  the 
dimension  of  the  parameter  space  Is  truncated.  The  resulting  parameter 
trade-offs  may  be  determined  by  an  inspection  of  the  parameter  resolution 
matrix. 


R  »  V  , 


(15) 


which  reduces  to  the  identity  matrix  in  the  case  of  perfect  resolution.  If 
during  test  Inversions,  the  parameter  resolution  is  deemed  unacceptable, 
adjustments  may  be  made  in  the  parameter  weights  or  the  singular  value 
cut-off,  and  the  inversion  may  be  repeated.  The  Importance  of  particular 
observations  or  residuals  in  the  inversion  may  be  determined  from  the 
information  density  matrix. 


This  may  indicate,  for  example,  that  the  parcimeter  changes  are  most  sensitive 
to  the  absolute  amplitude  residual  at  nearer  receivers,  and  to  the  absolute 
time  residual  at  more  distant  receivers.  The  information  density  matrix  from 
test  inversions  may  suggest  modifications  to  the  observation  weights. 

Finally,  since  most  of  the  partial  derivatives  in  A  are  determined 
numerically,  some  care  must  be  taken  in  determining  the  size  of  the  par^uneter 
perturbations  from  which  the  partials  are  computed.  In  addition,  particularly 
in  the  presence  of  local  minima  in  the  residual  space,  it  may  be  necessary  to 
constrain  or  damp  the  parameter  changes,  and  it  may  be  useful  to  utilize 
second  derivative  information.  Some  of  these  points  will  be  discussed  in  the 
examples  of  preliminary  inversions  to  follow. 

PRELIMINARY  ASSESSMENT  OF  STRUCTURE  INVERSION 

As  a  simple  initial  test  of  the  inversion  for  velocity  structure,  we 
generated  a  synthetic  data  set  for  an  explosion  source  in  a  halfspace.  The 
compressional  velocity  in  the  half space  was  6.0  km/s,  the  shear  velocity  was 
3.44  km/s  (a  Poisson  solid)  and  the  density  was  3.0  g/cm^.  The  source  was 
represented  by  a  Helraberger-Hadley  (1981)  RDP  with  K=7.0  s”\  B=1.0  and 
'^„,=  1 .0x10^®  cm^,  and  was  located  at  a  depth  of  1.0  km.  Vertical  component 
velocity  synthetics  were  computed  for  ranges  of  2,  4,  6,  8,  10  and  12  km. 

Since  the  object  of  the  test  inversion  is  the  halfspace  velocity,  we  constrain 


the  structure  to  be  a  single  layer  with  a  zero  velocity  gradient. 

In  such  a  simple  case,  it  is  instructive  to  investigate  the  residuals  as 
a  function  of  the  parauneter  values.  In  Figure  1  are  the  waveform  residuals, 
e^,  absolute  amplitude  residuals,  a^,  and  travel  time  residuals,  tj^,  plotted 
as  functions  of  receiver  distance  and  half space  velocity.  Velocities  from  5.0 
to  7.0  km/s  are  considered  in  increments  of  0.25  km/s.  Since  the  body 
waveforms  do  not  change  with  variations  in  the  halfspace  velocity,  the 
variations  in  waveform  residuals  are  due  to  the  truncation  of  the  synthetic 
waveform  as  it  is  shifted  with  respect  to  the  time  window.  The  waveform 
residual  is  zero  only  when  there  is  a  perfect  match:  at  the  correct  half space 
velocity  of  6.0  km/s.  Similarly,  the  absolute  ^unplitude  residual  does  not 
represent  the  difference  in  peak  amplitudes  between  observed  and  synthetic 
waveforms,  but  rather  the  difference  in  the  peak  values  of  their 
autocorrelations.  So,  systematic  amplitude  variations  due  to  velocity  changes 
are  superposed  with  truncation  effects  due  to  the  time  shifts  relative  to  the 
window.  Amplitude  residuals  may  take  on  positive  and  negative  values,  and  the 
sign  is  important  in  determining  the  partials.  Only  at  the  correct  velocity 
(6.0  km/s)  are  the  amplitude  residuals  at  all  of  the  receivers  zero.  The  most 
sensitive  measure  of  the  correct  halfspace  velocity  appears  to  be  the  travel 
time  residual.  This  residual  is  linearly  decreasing  at  all  receivers  for 
Increasing  halfspace  velocity  and,  of  course,  the  nvost  distant  receivers  are 
most  sensitive  to  variations  in  the  velocity.  The  travel  time  residual  is 
zero  at  all  receivers  for  the  correct  velocity. 

For  simple  problems  with  only  a  couple  of  free  parameters,  a 
grid-searching  approach  in  which  the  minimum  of  the  residuals  is  found  by 
stepping  through  the  paremeter  space,  would  be  quite  reasonable,  particularly 


Halfspace  Model  Residuals 
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Figure  1.  Waveform,  absolute  amplitude  and  travel  time  residuals  for  the 
halfspace  model  as  functions  of  receiver  distance  and  velocity.  The 
correct  velocity  is  6.0  km/s. 
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in  the  presence  of  local  minima.  This  is  essentially  the  procedure  used  to 
generate  Figure  1.  The  grid-searching  approach,  however,  becor.as  quite 
unwieldy  when  the  number  of  free  pareuneters  exceeds  three  or  four,  and  does 
not  readily  supply  Information  on  the  relative  resolution  of  the  parameters. 

In  anticipation  of  more  complex  problems,  therefore,  we  have  performed  the 
generalized  inverse  for  the  test  case  of  determining  the  half space  velocity 
from  a  starting  model  with  a  velocity  of  5.0  km/s.  If  we  were  to  invert  only 
the  travel-time  residuals.  Figure  1  indicates  that  the  Inversion  would  be 
linear.  However,  since  the  partlals  are  computed  by  a  linear  difference  of 
the  residuals,  the  curvature  in  the  waveform  and  eimplitude  residuals  requires 
an  Iterative  scheme.  A  10%  perturbation  in  velocity  is  used  to  compute  the 
partlals.  Accounting  for  the  dimensionality  of  the  residuals,  and  utilizing 
the  observation  that  travel  time  residuals  appear  to  be  most  sensitive  to 
halfspace  velocity,  we  define  the  covariance  matrix  of  the  observations  to  be 
dlagd.O,  1.0,  0.1)  for  the  waveform,  amplitude  and  time  residuals  at  each 
receiver,  respectively.  Since  there  is  only  one  parameter,  the  covariance  of 
the  velocity  is  set  to  1.0.  Five  iterations  were  allowed,  with  no  constraint 
or  damping  placed  on  the  pareuneter  changes.  The  results  are  tabulated  in 
Table  1,  and  the  waveforms  for  each  Iteration  are  plotted  in  Figure  2. 

After  only  two  iterations,  the  velocity  determined  is  within  0.5%  of  the 
correct  value,  and  convergence  is  complete  after  four  iterations.  The  quality 
of  fit  may  be  determined  by  inspecting  the  waveforms  in  Figure  2,  or  by  two 
objective  measures  listed  in  Table  1.  The  RMS  fit  is  defined  as  RMS  = 

/  r  /  3*m,  where  r  is  the  objective  function  (11)  which  is  the  sum  of  the 
squared  residuals  and  3*m  is  the  total  number  of  residuals.  The  least-squares 
error  is  LSE  =  1  A  Ap  -  Acj^,  and  is  a  measure  of  the  linearity  of  the 


Table  1  -  Halfspace  Velocity  Inversion  Result 


Iteration 

a  (km/s) 

RMS  fit 

Least-square 

Start 

5.000 

0.277 

1 

5.629 

0.137 

1.978 

2 

5.969 

0.130 

0.280 

3 

5.987 

0.056 

0.308 

4 

5.999 

0.001 

0.056 

5 

5.999 

0.012 

0.0002 

Data 

6.000 

inversion  at  each  iteration.  In  initial  inversions,  the  time  windows  were 
chosen  to  include  the  halfspace  Rayleigh  wave  at  each  receiver.  This, 
however,  causes  a  strong  local  minimum  to  be  formed  at  low  velocities  when  the 
"synthetic"  P  wave  correlates  with  the  "observed"  Rayleigh  wave.  Although  the 
Generalized  Ray  synthetics  are  exact  for  the  half space  Rayleigh  wave,  they  are 
generally  inadequate  for  modeling  surface  waves  for  more  complicated 
structures.  For  both  of  these  reasons,  the  inversion  time  windows  were 
shortened  to  exclude  the  Rayleigh  wave  where  it  is  separable  from  the  P  wave. 

For  a  somewhat  more  realistic  test,  and  in  order  to  test  the 
discretization  and  ray  summing  procedures,  a  synthetic  data  set  was  generated 
for  a  single  velocity  gradient.  The  corapresslonal  velocity  at  the  surface  was 
set  to  3.5  km/s  and  the  velocity  gradient  was  0.5  km/s/km.  The  shear  velocity 
was  constrained  to  that  of  a  Poisson  solid  (1//3  times  the  compressional 
velocity)  and  the  density  was  constrained  to  1/2  the  compressional  velocity. 
The  gradient  was  discretized  using  a  predominant  period  of  0.25  s,  which 
results  in  10  layers  in  the  top  15  km  with  layer  thickness  increasing  with 
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Figure  2.  Results  of  the  half space  velocity  Inversion.  Synthetic  seismograms  are  shown  for  the 

starting  mo<lel  an^l  five  iterations  at  each  receiver.  The  data  seismograms  (also  synthetic)  a 
plotted  below.  Arrows  above  the  data  waveforms  Indicate  the  inversion  time  windows  and  numbe 
to  the  right  are  the  peak  trace  amplitudes.  The  values  of  the  parameter  (halfspace  velocity) 
for  each  iteration,  and  the  correct  value  of  the  data,  are  shown  to  the  left. 
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depth.  All  upgoing  and  diving  direct  P  waves  and  the  surface-reflected  pP 
waves  were  considered,  but  culling  late  arriving  rays  resulted  in  18  rays 
being  summed  for  the  final  response.  Once  again,  synthetic  velocity 
seismograms  were  computed  for  ranges  of  2,  4,  6,  8,  10  and  12  km. 

Figures  3  and  4  illustrate  the  dependence  of  the  waveform,  absolute 
amplitude  and  travel  time  residuals  at  each  receiver  as  functions  of  the 
velocity  gradient  and  the  velocity  at  the  surface,  respectively.  As  shown  in 
Figure  3,  since  the  waveforms  at  the  nearer  receivers  are  dominated  by  the 
upgoing  P  wave,  which  is  relatively  insensitive  to  the  velocity  gradient, 
there  is  very  little  variation  in  their  residuals  over  the  ranges  of  gradient 
considered  (0.3  to  0.7  km/s/km).  On  the  other  hand,  the  waveforms  at  the 
farther  receivers  are  dominated  by  diving  turning  rays,  and  their  residuals 
(particularly  travel  time)  can  resolve  changes  in  the  gradient.  Overall,  it 
appears  that  velocity  gradient,  in  this  simple  example,  will  be  resolved  by 
travel  time  residuals  at  farther  receivers  and  eunplltude  residuals  at  all 
receivers.  The  residuals  due  to  variations  in  surface  velocity  (Figure  4)  are 
very  similar  to  those  for  the  halfspace  test  (Figure  1).  The  waveform 
residuals  have  a  broad  minimum,  the  travel  time  residuals  are  a  robust  measure 
of  the  correct  surface  velocity,  and  the  ^lmplltude  residuals  include 
variations  due  to  truncation  by  the  time  window.  In  comparing  Figures  3  and 
4,  it  appears  that  the  residuals  are  more  sensit  ve  to  variations  in  surface 
velocity  than  velocity  gradient.  We  will  have  to  take  this  into  account  in 
assigning  pareuneter  weights  in  the  inversion. 

Before  performing  an  inversion  for  surface  velocity  and  gradient,  it  is 
instructive  to  invert  for  gradient  only,  constraining  the  surface  velocity  to 
its  correct  value.  In  this  case,  we  define  the  covariance  matrix  of  the 


Gradient  Model  Residuals 


Gradient  (km/s^'Ton) 


Figure  3.  Waveform,  amplitude  and  time  residuals  for  the  gradient  model  as 
functions  of  receiver  distance  and  gradient.  The  correct  gradient  is  0 
km/s/km. 
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observations  according  to  the  dimensionality  of  the  residuals,  but  in 
addition,  we  wish  to  preferentially  minimize  the  waveform  residuals  since 
changes  in  gradient  appear  to  have  significant  effect  on  the  details  of  the 
waveform.  We  define  the  matrix  to  be  diag(0.01,  1.0,  0.1)  for  waveform, 
amplitude  and  time  residuals  at  each  receiver,  respectively.  Partial 
derivatives  are  computed  by  10%  perturbations  in  the  gradient,  beginning  with 
a  gradient  of  0.6  km/s/km.  The  correct  gradient  is  0.5  km/s/km.  Five 
iterations  were  allowed  with  no  damping  of  parameter  changes.  The  results  of 
the  Inversion  are  listed  in  Table  2  and  the  waveforms  for  each  iteration  are 
plotted  in  Figure  5. 

Convergence  to  the  solution  is  much  more  difficult  in  this  case.  An 
inspection  of  the  information  density  matrix  indicates  that  the  first 
iteration  is  primarily  controlled  by  the  waveform  and  travel  time  residuals  at 
the  12  km  receiver,  with  some  influence  by  the  waveform  and  amplitude 
residuals  at  6  km.  After  this  initial  step  toward  the  correct  solution,  the 
second  and  third  iterations  are  controlled  by  the  eunplitude  residuals  at  6  and 
10  km.  As  may  be  seen  in  Figure  3,  the  eunplitude  residuals  at  these  receivers 
vary  more  rapidly  for  gradients  slightly  less  than  0.5  km/s/km  than  for  those 
slightly  greater  than  0.5  km/s/km.  In  the  second  iteration,  the  partial 
derivatives  are  computed  based  on  the  difference  in  the  residuals  at  0.573  and 
0.469  km/3/)tm.  Even  though  the  latter  is  closer  to  the  correct  solution  (0.5 
km/s/km),  its  amplitude  residual  is  larger  in  absolute  value  than  the  former, 
so  the  partial  derivative  of  the  amplitude  residuals  causes  the  solution  to 
move  away  from  the  correct  gradient.  By  the  fourth  iteration,  the  solution 
has  moved  far  enough  away  that  the  lower  bracket  value  used  in  computing  the 
partials  is  approximately  0.5  km/s/km,  and  the  partial  is  unaffected  by  the 


Table  2  -  Gradient  Inversion  Results 


Iteration 

da 

Sz 

{ km/s/km) 

RMS  fit 

Least-square 

Start 

0.600 

0.102 

1 

0.522 

0.052 

1.015 

2 

0.529 

0.107 

0.125 

3 

0.560 

0.083 

2.128 

4 

0.497 

0.047 

0.080 

5 

0.471 

0.117 

0.088 

Data 

0.500 

changes  in  the  residuals  for  smaller  gradients.  It  is  clear  from  this 
discussion  that  some  information  on  the  second  derivatives  of  the  residuals  is 
needed  in  determining  the  parameter  changes. 

Finally,  we  have  performed  several  test  inversions  for  simultaneously 
determining  velocity  gradient  and  surface  velocity.  There  is  an  obvious 
trade-off  between  these  parameters,  indicating  the  need  to  truncate  small 
singular  values.  In  addition,  however,  the  parameters  should  be  weighted 
beyond  that  required  by  the  dimensionality  of  the  parameters.  Figures  3  and  4 
illustrated  that  the  Inversion  appears  to  be  more  sensitive  to  variations  in 
surface  velocity  than  gradient.  This  means  that  a  given  residual  value  may  be 
due  either  to  a  reasonable  variation  in  surface  velocity  or  a  relatively 
larger  variation  in  gradient.  Since  the  inversion  is  non-linear,  an  error  in 
surface  velocity  will  be  only  partially  corrected,  while  the  gradient  varies 
substantially  to  compensate.  Inversions  in  which  the  pareuneters  are  equally 
weighted,  or  in  which  gradient  is  preferentially  weighted  will  diverge.  We 
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Figure  5.  Results  of  the  inversion  for  velocity  gradient.  The  format  Is  the  same  as  in  Figure 
except  that  the  parameter  shown  on  the  left  is  velocity  gradient  in  km/s/km. 
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have  therefore  defined  the  covariance  matrix  of  the  parameters  such  that 
surface  velocity  is  preferentially  weighted  as  diagd.O,  0.1)  for  surface 
velocity  and  gradient,  respectively.  The  covariance  matrix  of  the 
observations  is  diag(0.1,  1.0,  0.1)  for  waveform,  amplitude  and  travel  time 
residuals,  respectively.  A  cut-off  was  established  so  that  one  of  the 
singular  values  could  be  truncated  to  avoid  large  variations  due  to  poor 
resolution.  Perturbations  of  10%  were  used  to  generate  the  partials.  Ten 
iterations  were  allowed  from  a  starting  model  with  surface  velocity  of  2.5 
)cm/s  and  a  gradient  of  0.4  km/a/km.  The  correct  solution  had  a  velocity  of 
3.5  km/s  and  a  gradient  of  0.5  km/s/km. 

The  results  of  the  inversion  are  listed  in  Table  3,  and  the  waveforms 
from  chosen  iterations  are  shown  in  Figure  6.  In  the  first  iteration,  both 
parameters  are  perfectly  resolved  and  the  change  is  toward  the  correct 
solution  for  both.  In  the  second  and  following  Iterations,  one  singular  value 
has  been  truncated.  The  weighting  scheme  has  ensured  that  the  remaining 
eigenvector  is  oriented  in  the  parameter  space  mostly  in  the  direction  of 
surface  velocity,  but  with  a  small  component  in  the  gradient.  This  means  that 
the  second  and  following  Iterations  consist  mainly  of  modifications  to  the 
surface  velocity,  with  only  small  (and  perhaps  incorrect)  modifications  to  the 
gradient.  This  is  also  indicated  by  the  diagonals  of  the  resolution  matrix, 
which  are  included  in  Table  3.  After  about  five  iterations,  the  model  is 
slightly  low  in  surface  velocity,  which  is  counteracted  by  a  slightly  high 
gradient.  The  waveforms  from  this  combination  are  close  enough  to  those  of 
the  correct  solution,  that  parameter  perturbations  are  relatively  small  in 
later  iterations.  Figure  6  shows  that  the  primary  effect  of  these  later 
perturbations  is  to  Improve  the  absolute  eunplitude,  since  travel  time  and 
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Iteration 

Cl  (km/s) 

* 

R 

(km/s Am) 

R 

RMS  fit 

LSE 

Start 

2.500 

0.400 

0.598 

1 

2.751 

1.000 

0.495 

1.000 

0.311 

25.9 

2 

2.994 

0.967 

0.509 

0.033 

0.238 

10.1 

3 

3.331 

0.921 

0.540 

0.079 

0.156 

1.01 

4 

3.453 

0.936 

0.550 

0.064 

0.139 

0.409 

5 

3.458 

0.842 

0.551 

0.158 

0.094 

0.369 

6 

3.437 

0.947 

0.549 

0.053 

0.144 

0.184 

7 

3.439 

0.947 

0.549 

0.053 

0.151 

0.406 

8 

3.432 

0.917 

0.549 

0.083 

0.140 

0.470 

9 

3.456 

0.967 

0.550 

0.033 

0.129 

0.363 

10 

3.466 

0.881 

0.551 

0.119 

0.078 

0.311 

Data 

3.500 

0.500 

* 


value  of  the  diagonal  of  the  resolution  matrix  for  this 


parameter. 


waveshape  are  already  fairly  well  fit.  This  indicates  that  it  may  be 
desirable  to  allow  the  covariance  matrix  of  the  observations  to  change  for 
later  iterations,  to  allow  the  amplitude  residual  more  influence  on  the 
solution. 

COMCLOSIONS 

We  have  introduced  a  method  for  the  simultaneous  Inversion  of  near-field 
waveforms  for  source  and  structure  pareuneters.  The  source  may  be  represented 
by  any  of  a  number  of  effective  source  functions.  The  structure  model  is 
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parameterized  as  a  layered  stack ,  in  which  the  free  par^u^^eters  are  the 
velocity  at  the  top  of  each  layer,  the  velocity  gradient  within  the  layer  and 
the  depth  to  the  top  of  the  layer.  For  each  observation,  the  residuals 
consist  of  the  normalized  cross-correlation  coefficient  between  the  data  and 
synthetics  (a  measure  of  waveform  fit),  the  difference  of  the  normalizations 
(a  measure  of  absolute  amplitude),  and  the  time  lag  to  the  peak  of  the 
cross-correlation  (a  measure  of  absolute  travel  time).  Numerical  partial 
derivatives  are  computed  using  asymptotic  Generalized  Ray  synthetics,  and  the 
problem  is  solved  using  an  iterative  generalized  Inverse. 

Test  Inversions  have  been  performed  for  simple  structure  models.  The 
Inversion  for  halfspace  velocity  is  quite  robust,  but  an  inversion  for  a 
single  gradient  Indicates  that  the  curvature  of  the  residuals  in  the  par^uneter 
space  causes  the  linear  difference  scheme  for  computing  the  partials  to  be 
Inadequate.  Information  on  the  second  derivatives  of  the  residuals  could 
solve  this  problem.  In  solving  for  both  surface  velocity  and  gradient,  these 
parameters  trade  off.  Without  adequate  weighting  of  the  pareuneters,  the 
inversion  diverges.  However,  utilizing  the  fact  that  the  inversion  is  more 
sensitive  to  surface  velocity  than  gradient,  and  truncating  singular  values  to 
damp  poorly  resolved  pareuneter  changes,  causes  the  inversion  to  slowly 
approach  the  correct  solution.  Preferential  weighting  of  the  different 
residuals  as  a  function  of  iteration  may  improve  convergence. 

Applying  the  simultaneous  inversion  for  source  and  structure  pareuneters 
to  actual  data  will  not  be  a  simple  task.  Noise  in  the  data  and  constraints 
in  the  parameter izatlons  will  compound  the  uncertainty  caused  by  pareuneter 
trade-offs.  By  adjusting  parameter  and  observation  weights,  and  progressing 
from  simple  to  more  complex  pariuneterizations,  it  should  be  possible  to  obtain 


a  solution  comparable  to  those  obtained  by  trial-and-error  modeling.  The 
inversion,  however,  provides  two  additional  benefits  that  trial-and-error 
modeling  cannot;  (1)  it  automatically  provides  a  quantitative  measure  of  the 
resolution  of  the  parameters  and  the  Importance  of  the  observations  in 
obtaining  the  solution,  and  (2)  it  avoids  the  tedious  and  ambiguous  process  of 
defining  a  detailed  structure  model  for  a  given  site  before  solving  for  source 
parameters.  Further  tests  will  investigate  the  trade-offs  in  multiple 
gradient  models,  as  well  as  the  effect  of  different  source  pareuneterizations 
on  the  structure  determination.  Following  this,  we  will  apply  the  inversion 
to  observed  near-field  data  from  Pahute  Mesa  and  Yucca  Flats,  which  we  have 
previously  modeled  using  trial-and-error  techniques.  If  this  proves  fruitful, 
the  next  step  will  be  to  invert  near-field  data  from  different  test  sites, 
such  as  the  GASBUGGY,  ROLISON  and  RIO  BLANCO  events. 
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